Implement Molodensky-Badekas transformation

[rouault]

This is a variation of Helmert 7 parameters where the rotation is not applied regarding the Earth centre, but with a reference point (in cartesian space). So X,Y,Z input coordinates are expressed regarding this Xp,Yp,Zp evaluation point, the Helmert rotation is applied, and then we re-add Xp,Yp,Zp, and then the translation terms of Helmert. So the existing PJ_helmert could just be enhanced with those Xp,Yp,Zp coordinates (one can convince oneself easily that when Xp=Yp=Zp=0, this is standard 7-parameter Helmert). Of course this is better documented in EPSG guidance 7.2

Used for example for https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:coordinateOperation:EPSG::1066

[kbevers]

The existing Helmert is already a complicated beast. My immediate thought is that it would be better to put this in it's own operation that leverages the code in PJ_helmert.c. Would that be possible without creating a big web of spaghetti?

[rouault]

My reasoning for including this in PJ_helmert is that it is should be just a matter of doing something like

diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index 7acbbb3..655e402 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -375,18 +375,18 @@ static XYZ helmert_forward_3d (LPZ lpz, PJ *P) {
 
     scale = 1 + Q->scale * 1e-6;
 
-    X = lpz.lam;
-    Y = lpz.phi;
-    Z = lpz.z;
+    X = lpz.lam - Q->refp.x;
+    Y = lpz.phi - Q->refp.y;
+    Z = lpz.z - Q->refp.z;
 
 
     point.xyz.x = scale * ( R00 * X  +   R01 * Y   +   R02 * Z);
     point.xyz.y = scale * ( R10 * X  +   R11 * Y   +   R12 * Z);
     point.xyz.z = scale * ( R20 * X  +   R21 * Y   +   R22 * Z);
 
-    point.xyz.x += Q->xyz.x;
-    point.xyz.y += Q->xyz.y;
-    point.xyz.z += Q->xyz.z;
+    point.xyz.x += Q->refp.x + Q->xyz.x;
+    point.xyz.y += Q->refp.y + Q->xyz.y;
+    point.xyz.z += Q->refp.z + Q->xyz.z;
 
     return point.xyz;
 }

(by the way the use of lpz.lam and lpz.phi here is super confusing as we are in the cartesian space in input space)

[kbevers]

My reasoning for including this in PJ_helmert is that it is should be just a matter of doing something like

I am aware that the math is simple this way. I am just worried about the growing number of parameters in Helmert.

The Molodensky-Badekas only uses the first 7 Helmert parameters, right? We could make a molobadekas (or whatever...), use the code you just proposed for both helmert and molobadekas, have the TRANSFORMATION(molobadekas, 0) setup only read the relevant parameters and use helmert_forward_3d() and helmert_reverse_3d()? Obviously a new function for reading the shared parameters would be useful in this approach.

I think this will also communicate better what transformation is performed.

(by the way the use of lpz.lam and lpz.phi here is super confusing as we are in the cartesian space in input space)

I agree. Feel free to fix it by adding a PJ_COORD for the input data.

[rouault]

@kbevers Your drafted plan sounds good to me. Indeed having a 18-parameter Helmert is not super user friendly :-)