[gdal-dev] Fwd: Transformation of image data between pixel/line and projected coordinates - using geolocation arrays, in both directions, without disk access
Even Rouault
even.rouault at spatialys.com
Tue May 16 12:34:49 PDT 2023
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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230516/1ef9e064/attachment-0001.htm>
More information about the gdal-dev
mailing list