<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
Hi!<br>
<br>
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:<br>
<ol>
<li>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 <a
href="https://lists.osgeo.org/pipermail/gdal-dev/2012-January/031502.html">here</a>
and a related test in the GDAL autotest suite (<a
href="https://github.com/OSGeo/gdal/blob/master/autotest/alg/transformgeoloc.py">here</a>).
However, I can´t get it to work for my specific case.</li>
<li>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.<br>
</li>
</ol>
<p>Here is the code I already have to convert a sample image from
cartesian to projected coordinates:</p>
<blockquote>
<p><font size="2">import os<br>
from tempfile import TemporaryDirectory<br>
from osgeo import gdal, osr<br>
import numpy as np<br>
from matplotlib import pyplot as plt<br>
<br>
<br>
# get some test data<br>
swath_data = np.random.randint(1, 100, (500, 400))<br>
lons, lats = np.meshgrid(np.linspace(3, 5, 500),<br>
np.linspace(40, 42, 400))<br>
<br>
with TemporaryDirectory() as td:<br>
p_lons_tmp = os.path.join(td, 'lons.tif')<br>
p_lats_tmp = os.path.join(td, 'lats.tif')<br>
p_data_tmp = os.path.join(td, 'data.tif')<br>
p_data_vrt = os.path.join(td, 'data.vrt')<br>
p_data_mapgeo_vrt = os.path.join(td, 'data_mapgeo.vrt')<br>
<br>
# save numpy arrays to temporary tif files<br>
for arr, path in zip((swath_data, lons, lats),
(p_data_tmp, p_lons_tmp, p_lats_tmp)):<br>
rows, cols = arr.shape<br>
drv = gdal.GetDriverByName('GTiff')<br>
ds = drv.Create(path, cols, rows, 1, gdal.GDT_Float64)<br>
ds.GetRasterBand(1).WriteArray(arr)<br>
del ds<br>
<br>
# add geolocation information to input data<br>
wgs84_wkt = osr.GetUserInputAsWKT('WGS84')<br>
utm_wkt = osr.GetUserInputAsWKT('EPSG:32632')<br>
ds = gdal.Translate(p_data_vrt, p_data_tmp, format='VRT')<br>
ds.SetMetadata(<br>
<br>
dict(<br>
SRS=wgs84_wkt,<br>
X_DATASET=p_lons_tmp,<br>
Y_DATASET=p_lats_tmp,<br>
X_BAND='1',<br>
Y_BAND='1',<br>
PIXEL_OFFSET='0',<br>
LINE_OFFSET='0',<br>
PIXEL_STEP='1',<br>
LINE_STEP='1'<br>
),<br>
'GEOLOCATION'<br>
)del ds<br>
<br>
# warp from geolocation arrays and read the result<br>
gdal.Warp(p_data_mapgeo_vrt, p_data_vrt, format='VRT',
geoloc=True,<br>
srcSRS=wgs84_wkt, dstSRS=utm_wkt)<br>
data_mapgeo = gdal.Open(p_data_mapgeo_vrt).ReadAsArray()<br>
<br>
# visualize input and output data<br>
fig, axes = plt.subplots(1, 4)<br>
for i, (arr, title) in enumerate(zip((swath_data, lons, lats,
data_mapgeo),<br>
('swath data', 'lons',
'lats', 'projected data'))):<br>
axes[i].imshow(arr, cmap='gray')<br>
axes[i].set_title(title)<br>
plt.tight_layout()<br>
plt.show()</font></p>
<p><font size="2"><br>
</font></p>
</blockquote>
<p>Any help would be highly appreciated!</p>
<p>Best,</p>
<p>Daniel Scheffler<br>
</p>
<br>
<pre class="moz-signature" cols="72">--
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: <a class="moz-txt-link-abbreviated moz-txt-link-freetext" href="mailto:daniel.scheffler@gfz-potsdam.de">daniel.scheffler@gfz-potsdam.de</a></pre>
</body>
</html>