<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
Thanks, Even, that works perfectly! Very nice to have a
GDAL-implementation for that!<br>
<br>
Best,<br>
Daniel<br>
<br>
<br>
<br>
<div class="moz-cite-prefix">Am 16.05.2023 um 21:34 schrieb Even
Rouault:<br>
</div>
<blockquote type="cite"
cite="mid:5ee98455-952d-cb55-e6c2-be7898804a00@spatialys.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<p>Daniel,</p>
<p>You rightly spotted <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/gdal/pull/6069"
moz-do-not-send="true">https://github.com/OSGeo/gdal/pull/6069</a>
as the enabler for that capability.<br>
</p>
<p>if your existing target dataset has a geolocation array
attached to do it, this should just be a matter of doing:<br>
</p>
<p>target_ds = gdal.Open( filename, gdal.GA_Update )<br>
</p>
<p>gdal.Warp(target_ds, source_ds, ... other options here ...)</p>
<p>If the target dataset doesn't have a geolocation array attached
to it, you can point to an external one with the
DST_GEOLOC_ARRAY tranformer option<br>
</p>
<p>gdal.Warp(target_ds, source_ds,
transformerOptions=["DST_METHOD=GEOLOC_ARRAY",
"DST_GEOLOC_ARRAY=/path/to/geoloc_dataset"], ... other options
here ...)</p>
<p>Even<br>
</p>
<div class="moz-cite-prefix">Le 16/05/2023 à 19:35, Daniel
Scheffler a écrit :<br>
</div>
<blockquote type="cite"
cite="mid:c3158573-ee48-2e4f-1516-ca80b4e6aa13@gfz-potsdam.de">
<meta http-equiv="Content-Type" content="text/html;
charset=UTF-8">
Hi!<br>
<br>
Some time ago, I asked for an inversion of gdal.Warp based on
GEOLOCATION arrays (longitude/latitude). Back then, my question
was: 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?<br>
<br>
In 11/2021, this was unfortunately not implemented yet. However,
is seems like in the meantime someting like this has been added:<br>
- <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/gdal/pull/5520"
moz-do-not-send="true">https://github.com/OSGeo/gdal/pull/5520</a><br>
- <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/gdal/pull/6069"
moz-do-not-send="true">https://github.com/OSGeo/gdal/pull/6069</a><br>
- <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/gdal/blob/c92b22d02c99eae0152f49595947fb3747ddc280/autotest/gcore/geoloc.py#L396"
moz-do-not-send="true">https://github.com/OSGeo/gdal/blob/c92b22d02c99eae0152f49595947fb3747ddc280/autotest/gcore/geoloc.py#L396</a><br>
<br>
But I am not quite sure if that is what I want. If so, how would
a Python implementation based on gdal.Warp look like? Is that
documented somewhere?<br>
<br>
Best,<br>
Daniel<br>
<br>
<br>
<br>
<div class="moz-cite-prefix">Am 30.11.2021 um 14:20 schrieb
Daniel Scheffler:<br>
</div>
<blockquote type="cite"
cite="mid:b1cf6b86-4922-ca9a-bdd4-ec4fff2883d7@gfz-potsdam.de">
<meta http-equiv="Content-Type" content="text/html;
charset=UTF-8">
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 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>
<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="moz-mime-attachment-header"></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>
<br>
<pre class="moz-signature" cols="72">--
Dr. Daniel Scheffler
Helmholtz Centre Potsdam
GFZ German Research Centre For Geosciences
Department 1 - Geodesy and Remote Sensing
Section 1.4 - Remote Sensing and Geoinformatics
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="moz-mime-attachment-header"></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">--
Dr. Daniel Scheffler
Helmholtz Centre Potsdam
GFZ German Research Centre For Geosciences
Department 1 - Geodesy and Remote Sensing
Section 1.4 - Remote Sensing and Geoinformatics
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>