[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
Tue May 16 10:35:40 PDT 2023
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230516/8b5724e7/attachment-0001.htm>
More information about the gdal-dev
mailing list