[gdal-dev] Custom vertical transformations

[Even Rouault]

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
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210820/bd444bba/attachment.html>