<html>
<head>
<meta http-equiv="Content-Type" content="text/html;
charset=windows-1252">
</head>
<body>
<p>Kristian,</p>
<p>your question kind of intersects
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/pipermail/gdal-dev/2021-August/054529.html">https://lists.osgeo.org/pipermail/gdal-dev/2021-August/054529.html</a></p>
<p>Two part answer:<br>
</p>
<p>A) I guess that there's an issue with axis order and units in
your pipeline. In particular if your DHYMSEA_1km_6150_542.tif has
a EPSG geographic CRS set, then the axis order of the coordinates
that will be provided to PROJ will be latitude, longitude in
degrees, so you likely need to change it to add an initial
conversion "+step +proj=unitconvert +xy_in=deg +xy_out=rad +step
+proj=axisswap +order=2,1"</p>
<p>And you likely also need to add a final symetric operation at the
end of the string: "+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1"</p>
<p>and adding -t_srs EPSG:XXXX will probably make GDAL happier too
(the -ct pipeline will be the one actually used even if you
specify -s_srs and -t_srs)</p>
<p>Don't forget to add -overwrite if you don't erase your target
file between each attempt.</p>
<p>B) Now... despite all the above, you won't get the result you
expect (you should get a no-op when that "works"). The reason is
that the adjustment of values using the vertical transform done in
GDAL is implemented as a kind of a hack. gdalwarp is fundamentally
a 2D coordinate transformation engine. The 3D mode actually works
by warping a PROJ vertical grid as an intermediate raster, and
adding the pixel values to the result of the 2D warping (hope that
makes sense). So that can't work for an arbitrary pipeline as
provided with low-level -ct, which is only used for the 2D
transformation of pixel coordinates. It would be certainly
desirable to have a way of having generic vertical transformation
of a grid but that won't be a 5-minute fix. Please file an
enhancement ticket about that, so that we can at least track that.
The time-dependent aspect should work once the previous bigger
issue is addressed.<br>
</p>
<p>Regarding your immediate issue, the workaround I see would be
that you use PROJ to create a regular grid that contains the
vertical correction, use gdalwarp to transform it to the exact
extent and resolution of the result of your 2D-only gdalwarp, and
finally use gdal_calc to sum the two.</p>
<p>Even<br>
</p>
<div class="moz-cite-prefix">Le 20/08/2021 à 16:37, Kristian Evers
via gdal-dev a écrit :<br>
</div>
<blockquote type="cite"
cite="mid:d75bf232b8f347c49b858cd484786632@sdfe.dk">
<meta http-equiv="Content-Type" content="text/html;
charset=windows-1252">
<meta name="Generator" content="Microsoft Word 15 (filtered
medium)">
<style>@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}a:link, span.MsoHyperlink
{mso-style-priority:99;
color:#0563C1;
text-decoration:underline;}a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:#954F72;
text-decoration:underline;}span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri",sans-serif;
color:windowtext;}.MsoChpDefault
{mso-style-type:export-only;}div.WordSection1
{page:WordSection1;}</style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
<div class="WordSection1">
<p class="MsoNormal"><span lang="DA">Hi list,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="DA"><o:p> </o:p></span></p>
<p class="MsoNormal">I am trying to apply a custom vertical
transformation to a grid using gdalwarp and am struggling a
bit to get the results I want. It’s entirely possible I am
trying something not within the capabilities of gdalwarp but I
am sure someone here can help me figure that out.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I am devising a deformation model
transformation that is meant to adjust terrain heights over a
certain time period. The general idea here is to predict the
changes in the terrain in the future based on various
information like isostatic uplift and land subsidence. I have
this information available as gridded ground motion
velocities. For now though I am just concerned about making a
proof-of-concept transformation and subsequent terrain
adjustment using PROJ and GDAL. I’ve got the PROJ part figured
out but need a bit of help to do what I want with GDAL.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">My transformation is a using the PROJ
“defmodel” operation. My proof-of-concept is available in this
repository:
<a href="https://github.com/kbevers/dk2100"
moz-do-not-send="true">https://github.com/kbevers/dk2100</a>.
For now it’s just a dummy that uses the same deformation grid
twice and applies a correction to the input height value. An
example looks like this:
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal" style="text-indent:36.0pt">echo 687071.4391
6210141.3267 0 2025 | cct +proj=pipeline +step +proj=utm
+zone=32 +inv +step \<o:p></o:p></p>
<p class="MsoNormal" style="text-indent:36.0pt"><span lang="DA">+proj=defmodel
+model=./defmodel.json +step +proj=utm +zone=32<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:36.0pt"><span lang="DA">
687071.4391 6210141.3267 6.0000 2025.0000<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="DA"><o:p> </o:p></span></p>
<p class="MsoNormal">The UTM steps surrounding the defmodel
accomodate that the grids I want to transform have horizontal
coordinates registered as EPSG:25832. Ideally I want this to
be a time-dependent transformation but for now I am just
working with simple constant offsets. <o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">When using the above transformation with
gdalwarp I am presented with an error:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">gdalwarp -ct
"+proj=pipeline +step +proj=utm +zone=32 +inv +step
+proj=defmodel +model=./defmodel.json +step +proj=utm
+zone=32" \ DHYMSEA_1km_6150_542.tif out.tif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">ERROR 1: Too
many points (2601 out of 2601) failed to transform, unable to
compute output bounds.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">The immediate question is “what am I
missing”? I suppose I either have some errors in the
transformation setup or am missing some command line switches.
Any guidance is appreciated. The next question is, if I turn
the above transformation into a time-depending transformation,
would I be able to use that with gdalwarp? I see GDAL 3.4 has
some switches for applying coordinate epochs but haven’t been
able to test them since I don’t have version 3.4 readily
available at the moment.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Best regards,<o:p></o:p></p>
<p class="MsoNormal">Kristian<o:p></o:p></p>
</div>
<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>