[gdal-dev] Fwd: Transformation of image data between pixel/line and projected coordinates - using geolocation arrays, in both directions, without disk access
Daniel Scheffler
daniel.scheffler at gfz-potsdam.de
Wed May 17 06:14:30 PDT 2023
Thanks, Even, that works perfectly! Very nice to have a
GDAL-implementation for that!
Best,
Daniel
Am 16.05.2023 um 21:34 schrieb Even Rouault:
>
> Daniel,
>
> You rightly spotted https://github.com/OSGeo/gdal/pull/6069 as the
> enabler for that capability.
>
> if your existing target dataset has a geolocation array attached to do
> it, this should just be a matter of doing:
>
> target_ds = gdal.Open( filename, gdal.GA_Update )
>
> gdal.Warp(target_ds, source_ds, ... other options here ...)
>
> If the target dataset doesn't have a geolocation array attached to it,
> you can point to an external one with the DST_GEOLOC_ARRAY tranformer
> option
>
> gdal.Warp(target_ds, source_ds,
> transformerOptions=["DST_METHOD=GEOLOC_ARRAY",
> "DST_GEOLOC_ARRAY=/path/to/geoloc_dataset"], ... other options here ...)
>
> Even
>
> Le 16/05/2023 à 19:35, Daniel Scheffler a écrit :
>> Hi!
>>
>> Some time ago, I asked for an inversion of gdal.Warp based on
>> GEOLOCATION arrays (longitude/latitude). Back then, my question was:
>> 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?
>>
>> In 11/2021, this was unfortunately not implemented yet. However, is
>> seems like in the meantime someting like this has been added:
>> - https://github.com/OSGeo/gdal/pull/5520
>> - https://github.com/OSGeo/gdal/pull/6069
>> -
>> https://github.com/OSGeo/gdal/blob/c92b22d02c99eae0152f49595947fb3747ddc280/autotest/gcore/geoloc.py#L396
>>
>> But I am not quite sure if that is what I want. If so, how would a
>> Python implementation based on gdal.Warp look like? Is that
>> documented somewhere?
>>
>> Best,
>> Daniel
>>
>>
>>
>> Am 30.11.2021 um 14:20 schrieb Daniel Scheffler:
>>> Ok thanks, too bad that this is not implemented. I think the
>>> inversion of this transformation would be a nice feature to be added
>>> in GDAL. It would be very useful to me (especially if it is
>>> accessible via the Python bindings) and would ease the
>>> implementation in a processing pipeline for the upcoming EnMAP
>>> hyperspectral satellite. Should I open a feature request in the GDAL
>>> issue tracker on GitHub?
>>>
>>>
>>> Am 30.11.2021 um 14:02 schrieb Even Rouault:
>>>>
>>>>
>>>> Le 30/11/2021 à 12:52, Daniel Scheffler a écrit :
>>>>> Thanks a lot for taking the time, Even, I got the transformation
>>>>> from cartesian to projected coordinates to work in memory with the
>>>>> GTiff driver. With MEM, NUMPY or VRT it does not work because
>>>>> these formats are either not readable from /vsimem/ or don´t have
>>>>> a regular file path which is needed to set the X_DATASET and
>>>>> Y_DATASET keys in the GEOLOCATION metadata.
>>>> Ah indeed for X_DATASET/Y_DATASET, you can't use a MEM or NUMPY
>>>> dataset. But a /vsimem/xxx dataset in another format should work.
>>>>>
>>>>> Regarding the inversion of this transformation, i.e., from
>>>>> projected coordinates to pixel/line:
>>>>>> The logic of GDALCreateGenImgProjTransformer2() around
>>>>>> https://github.com/OSGeo/gdal/blob/master/alg/gdaltransformer.cpp#L1825
>>>>>> which is for the source dataset should be ported a few lines
>>>>>> after for the target dataset. Probably with a new transformer
>>>>>> option to be able to point to an auxiliary dataset, such as a VRT
>>>>>> one of your example, to extract the geolocation metadata items
>>>>>> from it, that would be different from the target dataset itself,
>>>>>> because, except perhaps for the netCDF case, most GDAL datasets
>>>>>> that expose a GEOLOCATION metadata domain must be read-only.
>>>>>
>>>>> I don´t completely get what you mean here. To me, this sounds like
>>>>> there might be a way to do the inverted transformation using the
>>>>> C-API of GDAL.
>>>> No, I meant there's some missing code to do that.
>>>>> However, I am a Python developer and my C skills are a bit poor.
>>>>> Is there any way to use the Python bindings here?
>>>>>
>>>>> Kind regards,
>>>>> Daniel
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> Am 26.11.2021 um 12:38 schrieb Even Rouault:
>>>>>> Daniel,
>>>>>>>
>>>>>>> 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
>>>>>>> <https://lists.osgeo.org/pipermail/gdal-dev/2012-January/031502.html>
>>>>>>> and a related test in the GDAL autotest suite (here
>>>>>>> <https://github.com/OSGeo/gdal/blob/master/autotest/alg/transformgeoloc.py>).
>>>>>>> However, I can´t get it to work for my specific case.
>>>>>>>
>>>>>> There's no reason it won't work with a in-memory dataset.
>>>>>>
>>>>>> If you use a MEM dataset, then you need to provide the dataset
>>>>>> object itself as the input dataset of gdal.Warp() (the name of a
>>>>>> MEM dataset is completely ignored. a MEM dataset can't be opened,
>>>>>> just created).
>>>>>>
>>>>>> Similarly with a NUMPY dataset. With GTiff and /vsimem/, they
>>>>>> behave as a regular file. If you pass it by name as input of
>>>>>> gdal.Warp(), you need to make sure to close (ds = None typically)
>>>>>> the dataset before, so it is properly flushed and can be opened.
>>>>>> But you can also pass it as an object without that constraint.
>>>>>>
>>>>>> You can also create purely in-memory VRT files by assigning them
>>>>>> an empty name. Then of course you need to provide them as objects
>>>>>> to gdal.Warp()
>>>>>>
>>>>>>
>>>>>>> 1. 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.
>>>>>>>
>>>>>> The logic of GDALCreateGenImgProjTransformer2() around
>>>>>> https://github.com/OSGeo/gdal/blob/master/alg/gdaltransformer.cpp#L1825
>>>>>> which is for the source dataset should be ported a few lines
>>>>>> after for the target dataset. Probably with a new transformer
>>>>>> option to be able to point to an auxiliary dataset, such as a VRT
>>>>>> one of your example, to extract the geolocation metadata items
>>>>>> from it, that would be different from the target dataset itself,
>>>>>> because, except perhaps for the netCDF case, most GDAL datasets
>>>>>> that expose a GEOLOCATION metadata domain must be read-only.
>>>>>>>
>>>>>>> 1.
>>>>>>>
>>>>>>>
>>>>>>> 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.scheffler at gfz-potsdam.de
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> gdal-dev mailing list
>>>>>>> gdal-dev at lists.osgeo.org
>>>>>>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>>>> --
>>>>>> http://www.spatialys.com
>>>>>> My software is free, but my time generally not.
>>>>>
>>>>> --
>>>>>
>>>>> 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.scheffler at gfz-potsdam.de
>>>>>
>>>>> _______________________________________________
>>>>> gdal-dev mailing list
>>>>> gdal-dev at lists.osgeo.org
>>>>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>> --
>>>> http://www.spatialys.com
>>>> My software is free, but my time generally not.
>>>>
>>>> _______________________________________________
>>>> gdal-dev mailing list
>>>> gdal-dev at lists.osgeo.org
>>>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>
>>> --
>>>
>>> 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.scheffler at gfz-potsdam.de
>>>
>>> _______________________________________________
>>> gdal-dev mailing list
>>> gdal-dev at lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> --
>>
>> Dr. Daniel Scheffler
>>
>> Helmholtz Centre Potsdam
>> GFZ German Research Centre For Geosciences
>> Department 1 - Geodesy and Remote Sensing
>> Section 1.4 - Remote Sensing and Geoinformatics
>> Telegrafenberg, 14473 Potsdam, Germany
>>
>> Phone: +49 (0)331/288-1198
>> e-mail:daniel.scheffler at gfz-potsdam.de
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
--
Dr. Daniel Scheffler
Helmholtz Centre Potsdam
GFZ German Research Centre For Geosciences
Department 1 - Geodesy and Remote Sensing
Section 1.4 - Remote Sensing and Geoinformatics
Telegrafenberg, 14473 Potsdam, Germany
Phone: +49 (0)331/288-1198
e-mail:daniel.scheffler at gfz-potsdam.de
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230517/aa3b0dfe/attachment-0001.htm>
More information about the gdal-dev
mailing list