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.scheff...@gfz-potsdam.de
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
