Kristian,
your question kind of intersects
https://lists.osgeo.org/pipermail/gdal-dev/2021-August/054529.html
Two part answer:
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"
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"
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)
Don't forget to add -overwrite if you don't erase your target file
between each attempt.
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.
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.
Even
Le 20/08/2021 à 16:37, Kristian Evers via gdal-dev a écrit :
Hi list,
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.
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.
My transformation is a using the PROJ “defmodel” operation. My
proof-of-concept is available in this repository:
https://github.com/kbevers/dk2100 <https://github.com/kbevers/dk2100>.
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:
echo 687071.4391 6210141.3267 0 2025 | cct +proj=pipeline +step
+proj=utm +zone=32 +inv +step \
+proj=defmodel +model=./defmodel.json +step +proj=utm +zone=32
687071.4391 6210141.3267 6.0000 2025.0000
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.
When using the above transformation with gdalwarp I am presented with
an error:
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
ERROR 1: Too many points (2601 out of 2601) failed to transform,
unable to compute output bounds.
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.
Best regards,
Kristian
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
