<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
Daniel,<br>
<blockquote type="cite"
cite="mid:8e78eb30-6a89-104b-d2c4-f9ea8231aced@gfz-potsdam.de"> <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"
moz-do-not-send="true">here</a> and a related test in the
GDAL autotest suite (<a
href="https://github.com/OSGeo/gdal/blob/master/autotest/alg/transformgeoloc.py"
moz-do-not-send="true">here</a>). However, I can´t get it to
work for my specific case.</li>
</ol>
</blockquote>
<p>There's no reason it won't work with a in-memory dataset.</p>
<p>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).</p>
<p>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.</p>
<p>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()<br>
</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:8e78eb30-6a89-104b-d2c4-f9ea8231aced@gfz-potsdam.de">
<ol>
<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>
</blockquote>
The logic of <span class="pl-en">GDALCreateGenImgProjTransformer2()
around
<a class="moz-txt-link-freetext" href="https://github.com/OSGeo/gdal/blob/master/alg/gdaltransformer.cpp#L1825">https://github.com/OSGeo/gdal/blob/master/alg/gdaltransformer.cpp#L1825</a>
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.<br>
</span>
<blockquote type="cite"
cite="mid:8e78eb30-6a89-104b-d2c4-f9ea8231aced@gfz-potsdam.de">
<ol>
<li> </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" moz-do-not-send="true">daniel.scheffler@gfz-potsdam.de</a></pre>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<pre class="moz-quote-pre" wrap="">_______________________________________________
gdal-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/gdal-dev">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="http://www.spatialys.com">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</body>
</html>