[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 Nov 30 05:20:05 PST 2021
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20211130/b9aa3adb/attachment-0001.html>
More information about the gdal-dev
mailing list