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

[Daniel Scheffler]

Hi!

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.
 2. 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.

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20211125/7430bc1f/attachment.html>