[PROJ] question regarding vertical shift

[Even Rouault]

Mathieu,

> I’m hoping the mailing list can enlighten me regarding handling of vertical
> shift in PROJ.
> 
> Using this Switzerland vertical shift grid (
> https://www.dropbox.com/s/6uzpg3gr47q1p1k/CHGeo04.gtx?dl=0), if I type in
> the following command to adjust vertical height (100 meters here):
> 
> echo 7.9995 46.7949 100 | cct +proj=vgridshift +grids=CHGeo04.gtx
> 
> PROJ returns the following result:
> 
> 7.9995000000   46.7949000000       49.8406           inf
> 
> The adjusted height is 49.8 meters. This puzzles me as the grid’s value at
> the provided longitude and latitude (i.e.  7.9995 46.7949) is 50.179. I
> would have assumed that this meant a vertical adjustment of +50.179, not
> -50.179 as observed by the PROJ returned values.
> 
> Is PROJ’s result the correct one based on the grid? Or should an inverse
> transformation (i.e. cct -t) be used here? If so, what’s the rationale of
> whether to pick forward or inverse transformation?

The maths used by PROJ are explained here:
https://proj.org/operations/transformations/vgridshift.html#cmdoption-arg-multiplier

So by default PROJ will subtract the value from the grid to the input Z value. The rationale
is mostly historic and related to the fact +proj=vgridshift was used to emulate PROJ.4 style
+geoidgrids= parameter, which when applied for a target CRS in the forward direction
must go from ellipsoidal height to orthometric height, hence subtracting the value in
a geoid grid (we could also justify this default multiplier=-1 as in EPSG geoid-based
transformations are documented to go from the Geographic 3D CRS to the Vertical CRS, which
also means that the grid value must be subtracted)

So if you want to add the value of a vgridshift / go from orthometric height to ellipsoidal
height, add +multiplier=1  (or using the -I flag, but I find +multiplier=1 less confusing
as you directly specify the math)

vgridshift can also be used to do transformations between 2 vertical CRS, and if the
grid is documented to transfrom from A to B (by adding value of the grid) and you want
to transform from A to B, then you must also specify +multiplier=1

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com