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