[gdal-dev] Custom vertical transformations

[Kristian Evers]

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. 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210820/e02e1d31/attachment.html>