[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 03:52:17 PST 2021
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.
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. 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20211130/745c22eb/attachment-0001.html>
More information about the gdal-dev
mailing list