<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
Ok thanks, too bad that this is not implemented. I think the
inversion of this transformation would be a nice feature to be added
in GDAL. It would be very useful to me (especially if it is
accessible via the Python bindings) and would ease the
implementation in a processing pipeline for the upcoming EnMAP
hyperspectral satellite. Should I open a feature request in the GDAL
issue tracker on GitHub?<br>
<br>
<br>
<div class="moz-cite-prefix">Am 30.11.2021 um 14:02 schrieb Even
Rouault:<br>
</div>
<blockquote type="cite"
cite="mid:1791a553-3e0e-f5e6-659d-11ac6f708ab3@spatialys.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<p><br>
</p>
<div class="moz-cite-prefix">Le 30/11/2021 à 12:52, Daniel
Scheffler a écrit :<br>
</div>
<blockquote type="cite"
cite="mid:c2cf3a49-200b-91a6-dd81-1b3059859558@gfz-potsdam.de">
<meta http-equiv="content-type" content="text/html;
charset=UTF-8">
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>
</blockquote>
Ah indeed for X_DATASET/Y_DATASET, you can't use a MEM or NUMPY
dataset. But a /vsimem/xxx dataset in another format should work.<br>
<blockquote type="cite"
cite="mid:c2cf3a49-200b-91a6-dd81-1b3059859558@gfz-potsdam.de">
<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. </div>
</blockquote>
No, I meant there's some missing code to do that.<br>
<blockquote type="cite"
cite="mid:c2cf3a49-200b-91a6-dd81-1b3059859558@gfz-potsdam.de">
<div class="moz-forward-container">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>
<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>
<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>
<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" href="mailto:daniel.scheffler@gfz-potsdam.de">daniel.scheffler@gfz-potsdam.de</a></pre>
</body>
</html>