[gdal-dev] Transformation of image data between pixel/line and projected coordinates - using geolocation arrays, in both directions, without disk access

2021-11-25 Thread Daniel Scheffler 

Hi!

I am trying to convert image data from cartesian/image coordinates to 
projected coordinates AND vice versa using geolocation arrays in GDAL. I 
have two questions:


1. Since this transformation is part of a processing chain implemented
   in Python, I try to transform the data directly in-memory, i.e,
   without any disk access. This saves IO time and avoids permission
   errors when trying to write temporary data on Windows. How can this
   be done? I got correct results with the code below, however, only
   when I temporarily write the data to disk. I tried to write the data
   to /vsimem/ using the MEM, GTiff and NUMPY drivers. However,
   gdal.Warp can´t find the data there (FileNotFoundError). I think,
   also the gdal.Transformer class might be useful and I found an
   interesting thread on that here
   
   and a related test in the GDAL autotest suite (here
   ).
   However, I can´t get it to work for my specific case.
2. My second question is how I can invert the transformation, i.e., how
   can I transform an image with projected coordinates back to
   cartesian/image coordinates, given that a geolocation array tells
   GDAL where to put which pixel in the output? Background is a
   processing pipeline for satellite data where some processing steps
   are running in sensor geometry (image data as acquired by the
   sensor, without any geocoding and projection) and I need to provide
   corresponding AUX data which originally come with projected coordinates.

Here is the code I already have to convert a sample image from cartesian 
to projected coordinates:


   import os
   from tempfile import TemporaryDirectory
   from osgeo import gdal, osr
   import numpy as np
   from matplotlib import pyplot as plt


   # get some test data
   swath_data = np.random.randint(1, 100, (500, 400))
   lons, lats = np.meshgrid(np.linspace(3, 5, 500),
 np.linspace(40, 42, 400))

   with TemporaryDirectory() as td:
    p_lons_tmp = os.path.join(td, 'lons.tif')
    p_lats_tmp = os.path.join(td, 'lats.tif')
    p_data_tmp = os.path.join(td, 'data.tif')
    p_data_vrt = os.path.join(td, 'data.vrt')
    p_data_mapgeo_vrt = os.path.join(td, 'data_mapgeo.vrt')

    # save numpy arrays to temporary tif files
    for arr, path in zip((swath_data, lons, lats), (p_data_tmp,
   p_lons_tmp, p_lats_tmp)):
    rows, cols = arr.shape
    drv = gdal.GetDriverByName('GTiff')
    ds = drv.Create(path, cols, rows, 1, gdal.GDT_Float64)
    ds.GetRasterBand(1).WriteArray(arr)
    del ds

    # add geolocation information to input data
    wgs84_wkt = osr.GetUserInputAsWKT('WGS84')
    utm_wkt = osr.GetUserInputAsWKT('EPSG:32632')
    ds = gdal.Translate(p_data_vrt, p_data_tmp, format='VRT')
    ds.SetMetadata(

    dict(
    SRS=wgs84_wkt,
    X_DATASET=p_lons_tmp,
    Y_DATASET=p_lats_tmp,
    X_BAND='1',
    Y_BAND='1',
    PIXEL_OFFSET='0',
    LINE_OFFSET='0',
    PIXEL_STEP='1',
    LINE_STEP='1'
    ),
    'GEOLOCATION'
    )del ds

    # warp from geolocation arrays and read the result
    gdal.Warp(p_data_mapgeo_vrt, p_data_vrt, format='VRT', geoloc=True,
  srcSRS=wgs84_wkt, dstSRS=utm_wkt)
    data_mapgeo = gdal.Open(p_data_mapgeo_vrt).ReadAsArray()

   # visualize input and output data
   fig, axes = plt.subplots(1, 4)
   for i, (arr, title) in enumerate(zip((swath_data, lons, lats,
   data_mapgeo),
  ('swath data', 'lons', 'lats',
   'projected data'))):
    axes[i].imshow(arr, cmap='gray')
    axes[i].set_title(title)
   plt.tight_layout()
   plt.show()


Any help would be highly appreciated!

Best,

Daniel Scheffler


--

M.Sc. Geogr. Daniel Scheffler

Helmholtz Centre Potsdam
GFZ German Research Centre For Geosciences
Department 1 - Geodesy and Remote Sensing
Section 1.4  - Remote Sensing
Telegrafenberg, 14473 Potsdam, Germany

Phone:  +49 (0)331/288-1198
e-mail:daniel.scheff...@gfz-potsdam.de
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL COG - Problem with Overviews

2021-11-25 Thread Rahkonen Jukka (MML) 
Hi,

Please use “reply to all” for including the list as a recipient.
More details would lead to less guessing. Please add gdalinfo about the source 
and result (georeferencing is not important), tell how did you get the image of 
the original (which bands were selected, was any color enhancement applied). 
QGIS and a look at Layer properties – Symbology is good for that. Test with 
gdal_translate before going to Python and send us the command if you face the 
same trouble. If you cannot solve the problem yourself send a link to test 
image. It does not need to be your real image, any technically similar image is 
good and usually even better if you cut it smaller. A 200x200 pixel sized piece 
of the original image without georeferencing is one alternative that probably 
would not reveal secrets.

-Jukka Rahkonen-

Lähettäjä: Nicole Kamp 
Lähetetty: torstai 25. marraskuuta 2021 15.22
Vastaanottaja: Rahkonen Jukka (MML) 
Aihe: Re: [gdal-dev] GDAL COG - Problem with Overviews

Hello!
Thanks. I tried it with 0,1,2 - it did not work.


On Thu, Nov 25, 2021 at 1:49 PM Rahkonen Jukka (MML) 
mailto:jukka.rahko...@maanmittauslaitos.fi>>
 wrote:
Hi,

Maybe the band list be zero based? Try [0,1,2].

-Jukka Rahkonen-

Lähettäjä: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> 
Puolesta Nicole Kamp
Lähetetty: torstai 25. marraskuuta 2021 14.44
Vastaanottaja: gdal-dev@lists.osgeo.org
Aihe: [gdal-dev] GDAL COG - Problem with Overviews

Dear GDAL Developers,
I have an overview problem with my Cloud Optimized GeoTIFF.
I tried to translate a 16 bit, 4 bands ortho image to a COG 8bit, 3 bands image 
with JPEG compression.

That is the original image.
That is the result.

It works fine with LZW compression, but I need the JPEG compression. What am I 
doing wrong?

Thanks,
Niki


Here is my code.
for raster in glob.glob(str(parameters[self.INPUT_Pfad])+'/*.tif'):
fileInfo = QFileInfo(raster)
baseName = fileInfo.baseName()
file_name = str(parameters[self.INPUT_Pfad])+'/'+baseName+'.tif'
feedback.pushInfo(file_name)
output_file1 = str(parameters[self.OUTPUT_Pfad])+'/'+baseName+'.tif'
output_file2 = 
str(parameters[self.OUTPUT_Pfad])+'/'+baseName+'_cog.tif'

# Create COG
x_size = 6250
y_size = 5000
num_bands = 4

src_ds = gdal.Open(file_name)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Translate
translate_options = gdal.TranslateOptions(format = 'COG', 
outputType = gdal.GDT_Byte, bandList =  
   [1,2,3],  scaleParams=[''],
  creationOptions = 
['TILED=YES','COMPRESS=JPEG','BIGTIFF=IF_NEEDED',   
   'BLOCKXSIZE=512', 'BLOCKYSIZE=512', 
'RESAMPLING=NEAREST',   
 'OVERVIEWS=IGNORE_EXISTING', 'INTERLEAVE=PIXEL']
  )
gdal.Translate(output_file1, src_ds, options=translate_options)

translate_options = None
src_ds = None
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL COG - Problem with Overviews

2021-11-25 Thread Rahkonen Jukka (MML) 
Hi,

Maybe the band list be zero based? Try [0,1,2].

-Jukka Rahkonen-

Lähettäjä: gdal-dev  Puolesta Nicole Kamp
Lähetetty: torstai 25. marraskuuta 2021 14.44
Vastaanottaja: gdal-dev@lists.osgeo.org
Aihe: [gdal-dev] GDAL COG - Problem with Overviews

Dear GDAL Developers,
I have an overview problem with my Cloud Optimized GeoTIFF.
I tried to translate a 16 bit, 4 bands ortho image to a COG 8bit, 3 bands image 
with JPEG compression.

That is the original image.
That is the result.

It works fine with LZW compression, but I need the JPEG compression. What am I 
doing wrong?

Thanks,
Niki


Here is my code.
for raster in glob.glob(str(parameters[self.INPUT_Pfad])+'/*.tif'):
fileInfo = QFileInfo(raster)
baseName = fileInfo.baseName()
file_name = str(parameters[self.INPUT_Pfad])+'/'+baseName+'.tif'
feedback.pushInfo(file_name)
output_file1 = str(parameters[self.OUTPUT_Pfad])+'/'+baseName+'.tif'
output_file2 = 
str(parameters[self.OUTPUT_Pfad])+'/'+baseName+'_cog.tif'

# Create COG
x_size = 6250
y_size = 5000
num_bands = 4

src_ds = gdal.Open(file_name)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Translate
translate_options = gdal.TranslateOptions(format = 'COG', 
outputType = gdal.GDT_Byte, bandList =  
   [1,2,3],  scaleParams=[''],
  creationOptions = 
['TILED=YES','COMPRESS=JPEG','BIGTIFF=IF_NEEDED',   
   'BLOCKXSIZE=512', 'BLOCKYSIZE=512', 
'RESAMPLING=NEAREST',   
 'OVERVIEWS=IGNORE_EXISTING', 'INTERLEAVE=PIXEL']
  )
gdal.Translate(output_file1, src_ds, options=translate_options)

translate_options = None
src_ds = None
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL COG - Problem with Overviews

2021-11-25 Thread Nicole Kamp 
Dear GDAL Developers,
I have an overview problem with my Cloud Optimized GeoTIFF.
I tried to translate a 16 bit, 4 bands ortho image to a COG 8bit, 3 bands
image with JPEG compression.

That is the original image. 
That is the result. 

It works fine with LZW compression, but I need the JPEG compression. What
am I doing wrong?

Thanks,
Niki


Here is my code.
for raster in glob.glob(str(parameters[self.INPUT_Pfad])+'/*.tif'):
fileInfo = QFileInfo(raster)
baseName = fileInfo.baseName()
file_name = str(parameters[self.INPUT_Pfad])+'/'+baseName+'.tif'
feedback.pushInfo(file_name)
output_file1 =
str(parameters[self.OUTPUT_Pfad])+'/'+baseName+'.tif'
output_file2 =
str(parameters[self.OUTPUT_Pfad])+'/'+baseName+'_cog.tif'

# Create COG
x_size = 6250
y_size = 5000
num_bands = 4

src_ds = gdal.Open(file_name)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Translate
translate_options = gdal.TranslateOptions(format = 'COG',
outputType = gdal.GDT_Byte, bandList =
   [1,2,3],  scaleParams=[''],
  creationOptions =
['TILED=YES','COMPRESS=JPEG','BIGTIFF=IF_NEEDED',
'BLOCKXSIZE=512', 'BLOCKYSIZE=512',
'RESAMPLING=NEAREST',
  'OVERVIEWS=IGNORE_EXISTING', 'INTERLEAVE=PIXEL']
  )
gdal.Translate(output_file1, src_ds, options=translate_options)

translate_options = None
src_ds = None
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev