[PROJ] API migration - Offset on z-axis

[Even Rouault]

Kilian,

> Our team upgraded proj API from 4 to 6.
> 
> We noticed difference between the result when using projection in
> Switzerland (for instance EPSG:2056).
> 
> Here, I try to recap the difference between our implementation:
> 
> *Initial implementation (PROJ 4):*
> 
> // Input: x,y,z be coordinate expressed in EPSG:4978
> 
> projPJ pj_in;
> projPJ pj_out;
> 
> pj_in = pj_init_plus("+init=epsg:4978");
> pj_out = pj_init_plus("+init=epsg:2056");
> 
> pj_transform(pj_in, pj_out, 1, 1, x, y, z);
> 
> *Current implementation (PROJ 6)*
> 
> // Input: PJ_COORD coord expressed in EPSG:4978
> PJ_CONTEXT* context = projContextCreate();
> PJ* p1 = proj_create_crs_to_crs(context, "EPSG:4978", "EPSG:2056", NULL);
> PJ* proj = proj_normalize_for_visualization(context, p1);
> coord = proj_trans(proj , PJ_FWD, coord);
> 
> ---
> 
> 
> We noticed that the output using our current implementation yield the
> wrong result: there is an offset along the z-axis.
> 
> Are we doing something wrong?

No, you aren't doing something wrong. This is more a feature/bug of PROJ, and some 
ambiguity in how coordinate operations are presented in the EPSG database

First, what can help when you're investigating issue is to use "projinfo" to look at the 
pipelines. 

For example with latest PROJ master:

projinfo -s EPSG:4978 -t EPSG:2056  -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Conversion from WGS 84 (geocentric) to WGS 84 (geog2D) + Inverse of CH1903+ 
to WGS 84 (1) + Swiss Oblique Mercator 1995, 1 m, Europe - Liechtenstein and Switzerland

PROJ string:
+proj=pipeline +step +inv +proj=cart +ellps=WGS84 +step +proj=push +v_3 +step +proj=cart 
+ellps=WGS84 +step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 +step +inv 
+proj=cart +ellps=bessel +step +proj=pop +v_3 +step +proj=somerc 
+lat_0=46.9524055555556 +lon_0=7.43958333333333 +k_0=1 +x_0=2600000 +y_0=1200000 
+ellps=bessel

You can see that there are a push v_3 and pop v_3 steps surrounding the cartesian 
conversions and the Helmert transformation. Their effect is to avoid the ellipsoidal height to 
be modified during the Helmert transformation. 
This behaviour is based on the EPSG dataset presenting the Helmert transformation to be in 
the geographic 2D domain, due to the source & target CRS for which the "CH1903+ to WGS 
84 (1)" operation is presented being 2D themselves. So for an operation on 2D CRS, changing 
ellipsoidal height doesn't make sense as there's no ellipsoidal height for such CRS.
Later discussions with experts at EPSG revealed that the intent is however that when using 
the 3D CRS based on those 2D CRS, ellipsoidal height transformation should still be done.

So let's look when using EPSG:4979 (WGS 84 3D geographic) instead of EPSG:4978:

$ projinfo -s EPSG:4979 -t EPSG:2056  -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of CH1903+ to WGS 84 (1) + Swiss Oblique Mercator 1995, 1 m, Europe - 
Liechtenstein and Switzerland

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg 
+z_in=m +xy_out=rad +z_out=m +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 
+step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 +step +inv +proj=cart 
+ellps=bessel +step +proj=pop +v_3 +step +proj=somerc +lat_0=46.9524055555556 
+lon_0=7.43958333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel

Same behaviour as above due to the fact that even if EPSG:4979 is 3D, EPSG:2056 is 2D.
If adding the --3d switch to force its promotion to 3D:

$ projinfo -s EPSG:4979 -t EPSG:2056  -o PROJ --3d
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of CH1903+ to WGS 84 (1) + Swiss Oblique Mercator 1995, 1 m, Europe - 
Liechtenstein and Switzerland

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg 
+z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=WGS84 +step +proj=helmert 
+x=-674.374 +y=-15.056 +z=-405.346 +step +inv +proj=cart +ellps=bessel +step +proj=somerc 
+lat_0=46.9524055555556 +lon_0=7.43958333333333 +k_0=1 +x_0=2600000 +y_0=1200000 
+ellps=bessel

one can see that there is no more a push / pop v_3.

However if trying the --3d trick with EPSG:4978 this doesn't lead to the expected result 
currently. So that part of the behaviour is more a bug on PROJ side.

API wise, --3d behaviour is obtained with using proj_crs_promote_to_3D().

For your particular case, you can either use manually the corrected pipeline without the push 
/ pop v_3, or use a more 2-step approach with a first conversion between EPSG:4978 and 
EPSG:4979 chained with a transformation between EPSG:4979 and EPSG:2056 promoted to 
3D. The correct behaviour on ellipsoidal height when using a Geographic 3D CRS was fixed in 
"recent" PROJ 6 or 7 versions, so you'd better using the latest PROJ 7.1.0
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20200716/73afe5be/attachment-0001.html>