<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
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 <font size="2">X_DATASET
and </font><font size="2">Y_DATASET keys in the GEOLOCATION
metadata</font>.<br>
<div class="moz-forward-container"> <br>
<font size="2">Regarding the inversion of this transformation,
i.e., </font><font size="2">from projected coordinates to
pixel/line:</font><br>
<blockquote type="cite">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"
moz-do-not-send="true">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.</span></blockquote>
<br>
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?<br>
<br>
Kind regards,<br>
Daniel<br>
<br>
<br>
<br>
<br>
<div class="moz-cite-prefix">Am 26.11.2021 um 12:38 schrieb Even
Rouault:<br>
</div>
<blockquote type="cite"
cite="mid:2c9d72a2-3950-f846-b3e4-ee80500179a1@spatialys.com">
<meta http-equiv="Content-Type" content="text/html;
charset=UTF-8">
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"
moz-do-not-send="true">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> <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" 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 moz-txt-link-freetext" href="mailto:gdal-dev@lists.osgeo.org" moz-do-not-send="true">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" moz-do-not-send="true">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" moz-do-not-send="true">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</blockquote>
<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>
</div>
</body>
</html>