[gdal-dev] Fwd: Transformation of image data between pixel/line and projected coordinates - using geolocation arrays, in both directions, without disk access

Tue, 30 Nov 2021 03:52:37 -0800

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.scheff...@gfz-potsdam.de

_______________________________________________
gdal-dev mailing list
gdal-dev@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.scheff...@gfz-potsdam.de

_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev






















					Reply via email to