diff --git a/docs/source/operations/transformations/helmert.rst b/docs/source/operations/transformations/helmert.rst index 85e140af36..73449b40c0 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= - Translation of the x-axis given in meters. + Translation of the y-axis given in meters. .. option:: +z= diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst index 45c4a80bb6..d4bbf4e0e6 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 0000000000..b7d638bfcc --- /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= + + Translation of the x-axis given in meters. + +.. option:: +y= + + Translation of the y-axis given in meters. + +.. option:: +z= + + Translation of the z-axis given in meters. + +.. option:: +s= + + Scale factor given in ppm. + +.. option:: +rx= + + X-axis rotation given arc seconds. + +.. option:: +ry= + + Y-axis rotation given in arc seconds. + +.. option:: +rz= + + Z-axis rotation given in arc seconds. + +.. option:: +px= + + Coordinate along the x-axis of the reference point given in meters. + +.. option:: +py= + + Coordinate along the y-axis of the reference point given in meters. + +.. option:: +pz= + + 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 7acbbb3657..757cf950f7 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 8c72dd1fc2..3592dcc44c 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 6ac0c70ca1..13b77b0a33 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -410,6 +410,30 @@ operation proj=helmert transpose 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 -------------------------------------------------------------------------------