Merge pull request #1161 from rouault/molobadekas

Implement Molodensky-Badekas transform (fixes #1160)

6 files changed, 345 insertions, 74 deletions

diff --git a/docs/source/operations/transformations/helmert.rst b/docs/source/operations/transformations/helmert.rst
index 85e140af..73449b40 100644
--- a/docs/source/operations/transformations/helmert.rst
+++ b/docs/source/operations/transformations/helmert.rst
@@ -7,7 +7,7 @@ Helmert transform
 .. versionadded:: 5.0.0
 
 The Helmert transformation changes coordinates from one reference frame to
-anoether by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
+another by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
 14-parameter kinematic counterparts.
 
 
@@ -25,7 +25,7 @@ anoether by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
 | **Output type** | Cartesian coordinates                                             |
 +-----------------+-------------------------------------------------------------------+
 
-The Helmert transform, in all it's various incarnations, is used to perform reference
+The Helmert transform, in all its various incarnations, is used to perform reference
 frame shifts. The transformation operates in cartesian space. It can be used to transform
 planar coordinates from one datum to another, transform 3D cartesian
 coordinates from one static reference frame to another or it can be used to do fully
@@ -117,7 +117,7 @@ Parameters
 
 .. option:: +y=<value>
 
-    Translation of the x-axis given in meters.
+    Translation of the y-axis given in meters.
 
 .. option:: +z=<value>
 
diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst
index 45c4a80b..d4bbf4e0 100644
--- a/docs/source/operations/transformations/index.rst
+++ b/docs/source/operations/transformations/index.rst
@@ -16,5 +16,6 @@ systems are based on different datums.
    helmert
    horner
    molodensky
+   molobadekas
    hgridshift
    vgridshift
diff --git a/docs/source/operations/transformations/molobadekas.rst b/docs/source/operations/transformations/molobadekas.rst
new file mode 100644
index 00000000..b7d638bf
--- /dev/null
+++ b/docs/source/operations/transformations/molobadekas.rst
@@ -0,0 +1,147 @@
+.. _molobadekas:
+
+================================================================================
+Molodensky-Badekas transform
+================================================================================
+
+.. versionadded:: 6.0.0
+
+The Molodensky-Badekas transformation changes coordinates from one reference frame to
+another by means of a 10-parameter shift.
+
+.. note::
+
+    It should not be confused with the :ref:`Molodensky` transform which
+    operates directly in the geodetic coordinates. Molodensky-Badekas can rather
+    be seen as a variation of :ref:`Helmert`
+
++-----------------+-------------------------------------------------------------------+
+| **Alias**       | molobadekas                                                       |
++-----------------+-------------------------------------------------------------------+
+| **Domain**      | 3D                                                                |
++-----------------+-------------------------------------------------------------------+
+| **Input type**  | Cartesian coordinates                                             |
++-----------------+-------------------------------------------------------------------+
+| **Output type** | Cartesian coordinates                                             |
++-----------------+-------------------------------------------------------------------+
+
+The Molodensky-Badekas transformation is a variation of the :ref:`Helmert` where
+the rotational terms are not directly applied to the ECEF coordinates, but on
+cartesian coordinates relative to a reference point (usually close to Earth surface,
+and to the area of use of the transformation). When ``px`` = ``py`` = ``pz`` = 0,
+this is equivalent to a 7-parameter Helmert transformation.
+
+Example
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+Transforming coordinates from La Canoa to REGVEN:
+
+::
+
+    proj=molobadekas convention=coordinate_frame
+           x=-270.933 y=115.599 z=-360.226 rx=-5.266 ry=-1.238 rz=2.381
+           s=-5.109 px=2464351.59 py=-5783466.61 pz=974809.81
+
+
+Parameters
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. note::
+
+    All parameters (except convention) are optional but at least one should be
+    used, otherwise the operation will return the coordinates unchanged.
+
+.. option:: +convention=coordinate_frame/position_vector
+
+    Indicates the convention to express the rotational terms when a 3D-Helmert /
+    7-parameter more transform is involved.
+
+    The two conventions are equally popular and a frequent source of confusion.
+    The coordinate frame convention is also described as an clockwise
+    rotation of the coordinate frame. It corresponds to EPSG method code
+    1034 (in the geocentric domain) or 9636 (in the geographic domain)
+    The position vector convention is also described as an anticlockwise
+    (counter-clockwise) rotation of the coordinate frame.
+    It corresponds to as EPSG method code 1061 (in the geocentric domain) or
+    1063 (in the geographic domain).
+
+    The result obtained with parameters specified in a given convention
+    can be obtained in the other convention by negating the rotational parameters
+    (``rx``, ``ry``, ``rz``)
+
+.. option:: +x=<value>
+
+    Translation of the x-axis given in meters.
+
+.. option:: +y=<value>
+
+    Translation of the y-axis given in meters.
+
+.. option:: +z=<value>
+
+    Translation of the z-axis given in meters.
+
+.. option:: +s=<value>
+
+    Scale factor given in ppm.
+
+.. option:: +rx=<value>
+
+    X-axis rotation given arc seconds.
+
+.. option:: +ry=<value>
+
+    Y-axis rotation given in arc seconds.
+
+.. option:: +rz=<value>
+
+    Z-axis rotation given in arc seconds.
+
+.. option:: +px=<value>
+
+    Coordinate along the x-axis of the reference point given in meters.
+
+.. option:: +py=<value>
+
+    Coordinate along the y-axis of the reference point given in meters.
+
+.. option:: +pz=<value>
+
+    Coordinate along the z-axis of the reference point given in meters.
+
+Mathematical description
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+In the *Position Vector* convention, we define :math:`R_x = radians \left( rx \right)`,
+:math:`R_z = radians \left( ry \right)` and :math:`R_z = radians \left( rz \right)`
+
+In the *Coordinate Frame* convention, :math:`R_x = - radians \left( rx \right)`,
+:math:`R_z = - radians \left( ry \right)` and :math:`R_z = - radians \left( rz \right)`
+
+.. math::
+    :label: 10param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^{output} =
+        \begin{bmatrix}
+            T_x + P_x\\
+            T_y + P_y\\
+            T_z + P_z\\
+        \end{bmatrix} +
+        \left(1 + s \times 10^{-6}\right)
+        \begin{bmatrix}
+             1  & -R_z  &  R_y \\
+             Rz &  1    & -R_x \\
+            -Ry &  R_x  &  1   \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            X^{input} - P_x\\
+            Y^{input} - P_y\\
+            Z^{input} - P_z\\
+        \end{bmatrix}
+    \end{align}
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index 7acbbb36..757cf950 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -10,6 +10,10 @@
   Implements 3(6)-, 4(8) and 7(14)-parameter Helmert transformations for
   3D data.
 
+  Also incorporates Molodensky-Badekas variant of 7-parameter Helmert
+  transformation, where the rotation is not applied regarding the centre
+  of the spheroid, but given a reference point.
+
   Primarily useful for implementation of datum shifts in transformation
   pipelines.
 
@@ -17,7 +21,8 @@
 
 Thomas Knudsen, thokn@sdfe.dk, 2016-05-24/06-05
 Kristian Evers, kreve@sdfe.dk, 2017-05-01
-Last update: 2017-05-15
+Even Rouault, even.roault@spatialys.com
+Last update: 2018-10-26
 
 ************************************************************************
 * Copyright (c) 2016, Thomas Knudsen / SDFE
@@ -52,6 +57,7 @@ Last update: 2017-05-15
 #include "geocent.h"
 
 PROJ_HEAD(helmert, "3(6)-, 4(8)- and 7(14)-parameter Helmert shift");
+PROJ_HEAD(molobadekas, "Molodensky-Badekas transformation");
 
 static XYZ helmert_forward_3d (LPZ lpz, PJ *P);
 static LPZ helmert_reverse_3d (XYZ xyz, PJ *P);
@@ -66,6 +72,7 @@ struct pj_opaque_helmert {
     XYZ xyz;
     XYZ xyz_0;
     XYZ dxyz;
+    XYZ refp;
     PJ_OPK opk;
     PJ_OPK opk_0;
     PJ_OPK dopk;
@@ -375,16 +382,16 @@ 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.x += Q->xyz.x; /* for Molodensky-Badekas, Q->xyz already incorporates the Q->refp offset */
     point.xyz.y += Q->xyz.y;
     point.xyz.z += Q->xyz.z;
 
@@ -421,9 +428,9 @@ static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) {
     Z = (xyz.z - Q->xyz.z) / scale;
 
     /* Inverse rotation through transpose multiplication */
-    point.xyz.x  =  ( R00 * X   +   R10 * Y   +   R20 * Z);
-    point.xyz.y  =  ( R01 * X   +   R11 * Y   +   R21 * Z);
-    point.xyz.z  =  ( R02 * X   +   R12 * Y   +   R22 * Z);
+    point.xyz.x  =  ( R00 * X   +   R10 * Y   +   R20 * Z) + Q->refp.x;
+    point.xyz.y  =  ( R01 * X   +   R11 * Y   +   R21 * Z) + Q->refp.y;
+    point.xyz.z  =  ( R02 * X   +   R12 * Y   +   R22 * Z) + Q->refp.z;
 
     return point.lpz;
 }
@@ -467,30 +474,108 @@ static PJ_COORD helmert_reverse_4d (PJ_COORD point, PJ *P) {
 /* Arcsecond to radians */
 #define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0)
 
-/***********************************************************************/
-PJ *TRANSFORMATION(helmert, 0) {
-/***********************************************************************/
+
+static PJ* init_helmert_six_parameters(PJ* P) {
     struct pj_opaque_helmert *Q = pj_calloc (1, sizeof (struct pj_opaque_helmert));
     if (0==Q)
         return pj_default_destructor (P, ENOMEM);
     P->opaque = (void *) Q;
 
-    P->fwd4d  = helmert_forward_4d;
-    P->inv4d  = helmert_reverse_4d;
-    P->fwd3d  = helmert_forward_3d;
-    P->inv3d  = helmert_reverse_3d;
-    P->fwd    = helmert_forward;
-    P->inv    = helmert_reverse;
-
     /* In most cases, we work on 3D cartesian coordinates */
     P->left  = PJ_IO_UNITS_CARTESIAN;
     P->right = PJ_IO_UNITS_CARTESIAN;
-    /* But in the 2D case, the coordinates are projected */
+
+    /* Translations */
+    if (pj_param (P->ctx, P->params, "tx").i)
+        Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
+
+    if (pj_param (P->ctx, P->params, "ty").i)
+        Q->xyz_0.y = pj_param (P->ctx, P->params, "dy").f;
+
+    if (pj_param (P->ctx, P->params, "tz").i)
+        Q->xyz_0.z = pj_param (P->ctx, P->params, "dz").f;
+
+    /* Rotations */
+    if (pj_param (P->ctx, P->params, "trx").i)
+        Q->opk_0.o = pj_param (P->ctx, P->params, "drx").f * ARCSEC_TO_RAD;
+
+    if (pj_param (P->ctx, P->params, "try").i)
+        Q->opk_0.p = pj_param (P->ctx, P->params, "dry").f * ARCSEC_TO_RAD;
+
+    if (pj_param (P->ctx, P->params, "trz").i)
+        Q->opk_0.k = pj_param (P->ctx, P->params, "drz").f * ARCSEC_TO_RAD;
+
+    /* Use small angle approximations? */
+    if (pj_param (P->ctx, P->params, "bexact").i)
+        Q->exact = 1;
+
+    return P;
+}
+
+
+static PJ* read_convention(PJ* P) {
+
+    struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *)P->opaque;
+
+    /* In case there are rotational terms, we require an explicit convention
+     * to be provided. */
+    if (!Q->no_rotation) {
+        const char* convention = pj_param (P->ctx, P->params, "sconvention").s;
+        if( !convention ) {
+            proj_log_error (P, "helmert: missing 'convention' argument");
+            return pj_default_destructor (P, PJD_ERR_MISSING_ARGS);
+        }
+        if( strcmp(convention, "position_vector") == 0 ) {
+            Q->is_position_vector = 1;
+        }
+        else if( strcmp(convention, "coordinate_frame") == 0 ) {
+            Q->is_position_vector = 0;
+        }
+        else {
+            proj_log_error (P, "helmert: invalid value for 'convention' argument");
+            return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
+        }
+
+        /* historically towgs84 in PROJ has always been using position_vector
+         * convention. Accepting coordinate_frame would be confusing. */
+        if (pj_param_exists (P->params, "towgs84")) {
+            if( !Q->is_position_vector ) {
+                proj_log_error (P, "helmert: towgs84 should only be used with "
+                                "convention=position_vector");
+                return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
+            }
+        }
+    }
+
+    return P;
+}
+
+
+/***********************************************************************/
+PJ *TRANSFORMATION(helmert, 0) {
+/***********************************************************************/
+
+    struct pj_opaque_helmert *Q;
+
+    if( !init_helmert_six_parameters(P) ) {
+        return 0;
+    }
+
+    /* In the 2D case, the coordinates are projected */
     if (pj_param_exists (P->params, "theta")) {
         P->left  = PJ_IO_UNITS_PROJECTED;
         P->right = PJ_IO_UNITS_PROJECTED;
     }
 
+    P->fwd4d  = helmert_forward_4d;
+    P->inv4d  = helmert_reverse_4d;
+    P->fwd3d  = helmert_forward_3d;
+    P->inv3d  = helmert_reverse_3d;
+    P->fwd    = helmert_forward;
+    P->inv    = helmert_reverse;
+
+    Q = (struct pj_opaque_helmert *)P->opaque;
+
     /* Detect obsolete transpose flag and error out if found */
     if (pj_param (P->ctx, P->params, "ttranspose").i) {
         proj_log_error (P, "helmert: 'transpose' argument is no longer valid. "
@@ -517,26 +602,6 @@ PJ *TRANSFORMATION(helmert, 0) {
             Q->scale_0 = (P->datum_params[6] - 1) * 1e6;
     }
 
-    /* Translations */
-    if (pj_param (P->ctx, P->params, "tx").i)
-        Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
-
-    if (pj_param (P->ctx, P->params, "ty").i)
-        Q->xyz_0.y = pj_param (P->ctx, P->params, "dy").f;
-
-    if (pj_param (P->ctx, P->params, "tz").i)
-        Q->xyz_0.z = pj_param (P->ctx, P->params, "dz").f;
-
-    /* Rotations */
-    if (pj_param (P->ctx, P->params, "trx").i)
-        Q->opk_0.o = pj_param (P->ctx, P->params, "drx").f * ARCSEC_TO_RAD;
-
-    if (pj_param (P->ctx, P->params, "try").i)
-        Q->opk_0.p = pj_param (P->ctx, P->params, "dry").f * ARCSEC_TO_RAD;
-
-    if (pj_param (P->ctx, P->params, "trz").i)
-        Q->opk_0.k = pj_param (P->ctx, P->params, "drz").f * ARCSEC_TO_RAD;
-
     if (pj_param (P->ctx, P->params, "ttheta").i) {
         Q->theta_0 = pj_param (P->ctx, P->params, "dtheta").f * ARCSEC_TO_RAD;
         Q->fourparam = 1;
@@ -585,10 +650,6 @@ PJ *TRANSFORMATION(helmert, 0) {
     if (pj_param(P->ctx, P->params, "tt_obs").i)
         Q->t_obs = pj_param (P->ctx, P->params, "dt_obs").f;
 
-    /* Use small angle approximations? */
-    if (pj_param (P->ctx, P->params, "bexact").i)
-        Q->exact = 1;
-
     Q->xyz    =  Q->xyz_0;
     Q->opk    =  Q->opk_0;
     Q->scale  =  Q->scale_0;
@@ -599,34 +660,8 @@ PJ *TRANSFORMATION(helmert, 0) {
         Q->no_rotation = 1;
     }
 
-    /* In case there are rotational terms, we require an explicit convention
-     * to be provided. */
-    if (!Q->no_rotation) {
-        const char* convention = pj_param (P->ctx, P->params, "sconvention").s;
-        if( !convention ) {
-            proj_log_error (P, "helmert: missing 'convention' argument");
-            return pj_default_destructor (P, PJD_ERR_MISSING_ARGS);
-        }
-        if( strcmp(convention, "position_vector") == 0 ) {
-            Q->is_position_vector = 1;
-        }
-        else if( strcmp(convention, "coordinate_frame") == 0 ) {
-            Q->is_position_vector = 0;
-        }
-        else {
-            proj_log_error (P, "helmert: invalid value for 'convention' argument");
-            return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
-        }
-
-        /* historically towgs84 in PROJ has always been using position_vector
-         * convention. Accepting coordinate_frame would be confusing. */
-        if (pj_param_exists (P->params, "towgs84")) {
-            if( !Q->is_position_vector ) {
-                proj_log_error (P, "helmert: towgs84 should only be used with "
-                                "convention=position_vector");
-                return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
-            }
-        }
+    if( !read_convention(P) ) {
+        return 0;
     }
 
     /* Let's help with debugging */
@@ -653,3 +688,66 @@ PJ *TRANSFORMATION(helmert, 0) {
 
     return P;
 }
+
+
+/***********************************************************************/
+PJ *TRANSFORMATION(molobadekas, 0) {
+/***********************************************************************/
+
+    struct pj_opaque_helmert *Q;
+
+    if( !init_helmert_six_parameters(P) ) {
+        return 0;
+    }
+
+    P->fwd3d  = helmert_forward_3d;
+    P->inv3d  = helmert_reverse_3d;
+
+    Q = (struct pj_opaque_helmert *)P->opaque;
+
+    /* Scale */
+    if (pj_param (P->ctx, P->params, "ts").i) {
+        Q->scale_0 = pj_param (P->ctx, P->params, "ds").f;
+    }
+
+    Q->opk    =  Q->opk_0;
+    Q->scale  =  Q->scale_0;
+
+    if( !read_convention(P) ) {
+        return 0;
+    }
+
+    /* Reference point */
+    if (pj_param (P->ctx, P->params, "tpx").i)
+        Q->refp.x = pj_param (P->ctx, P->params, "dpx").f;
+
+    if (pj_param (P->ctx, P->params, "tpy").i)
+        Q->refp.y = pj_param (P->ctx, P->params, "dpy").f;
+
+    if (pj_param (P->ctx, P->params, "tpz").i)
+        Q->refp.z = pj_param (P->ctx, P->params, "dpz").f;
+
+
+    /* Let's help with debugging */
+    if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
+        proj_log_debug(P, "Molodensky-Badekas parameters:");
+        proj_log_debug(P, "x=  %8.5f  y=  %8.5f  z=  %8.5f", Q->xyz_0.x, Q->xyz_0.y, Q->xyz_0.z);
+        proj_log_debug(P, "rx= %8.5f  ry= %8.5f  rz= %8.5f",
+                Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD);
+        proj_log_debug(P, "s=  %8.5f  exact=%d%s", Q->scale, Q->exact,
+                       Q->is_position_vector ? "  convention=position_vector" :
+                       "  convention=coordinate_frame");
+        proj_log_debug(P, "px= %8.5f  py= %8.5f  pz= %8.5f",   Q->refp.x, Q->refp.y, Q->refp.z);
+    }
+
+    /* as an optimization, we incorporate the refp in the translation terms */
+    Q->xyz_0.x +=  Q->refp.x;
+    Q->xyz_0.y +=  Q->refp.y;
+    Q->xyz_0.z +=  Q->refp.z;
+
+    Q->xyz    =  Q->xyz_0;
+
+    build_rot_matrix(P);
+
+    return P;
+}
diff --git a/src/pj_list.h b/src/pj_list.h
index 8c72dd1f..3592dcc4 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -94,6 +94,7 @@ PROJ_HEAD(mil_os, "Miller Oblated Stereographic")
 PROJ_HEAD(mill, "Miller Cylindrical")
 PROJ_HEAD(misrsom, "Space oblique for MISR")
 PROJ_HEAD(moll, "Mollweide")
+PROJ_HEAD(molobadekas, "Molodensky-Badekas transform") /* implemented in PJ_helmert.c */
 PROJ_HEAD(molodensky, "Molodensky transform")
 PROJ_HEAD(murd1, "Murdoch I")
 PROJ_HEAD(murd2, "Murdoch II")
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 6ac0c70c..13b77b0a 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -411,6 +411,30 @@ expect     failure   errno invalid_arg
 
 
 -------------------------------------------------------------------------------
+Molodensky-Badekas from IOGP Guidance 7.2, Transformation from La Canoa to REGVEN
+between geographic 2D coordinate reference systems (EPSG Dataset transformation code 1771).
+Here just taking the Cartesian step of the transformation.
+-------------------------------------------------------------------------------
+operation  proj=molobadekas convention=coordinate_frame
+           x=-270.933 y=115.599 z=-360.226 rx=-5.266 ry=-1.238 rz=2.381
+           s=-5.109 px=2464351.59 py=-5783466.61 pz=974809.81
+-------------------------------------------------------------------------------
+tolerance  1 cm
+roundtrip  1
+accept     2550408.96 -5749912.26 1054891.11
+expect     2550138.45 -5749799.87 1054530.82
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+Test error cases of molobadekas
+-------------------------------------------------------------------------------
+
+# Missing convention
+operation  proj=molobadekas
+expect     failure   errno missing_arg
+
+
+-------------------------------------------------------------------------------
 geocentric latitude
 -------------------------------------------------------------------------------
 operation  proj=geoc ellps=GRS80