From df574ae332d57f556fd56314883b3354cab1d0ff Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 19 Dec 2018 13:00:37 +0100 Subject: cpp conversion: remove useless pj_, PJ_ and proj_ filename prefixes --- src/transformations/PJ_affine.cpp | 250 ----------- src/transformations/PJ_deformation.cpp | 325 -------------- src/transformations/PJ_helmert.cpp | 755 --------------------------------- src/transformations/PJ_hgridshift.cpp | 133 ------ src/transformations/PJ_horner.cpp | 513 ---------------------- src/transformations/PJ_molodensky.cpp | 329 -------------- src/transformations/PJ_vgridshift.cpp | 144 ------- src/transformations/affine.cpp | 250 +++++++++++ src/transformations/deformation.cpp | 325 ++++++++++++++ src/transformations/helmert.cpp | 755 +++++++++++++++++++++++++++++++++ src/transformations/hgridshift.cpp | 133 ++++++ src/transformations/horner.cpp | 513 ++++++++++++++++++++++ src/transformations/molodensky.cpp | 329 ++++++++++++++ src/transformations/vgridshift.cpp | 144 +++++++ 14 files changed, 2449 insertions(+), 2449 deletions(-) delete mode 100644 src/transformations/PJ_affine.cpp delete mode 100644 src/transformations/PJ_deformation.cpp delete mode 100644 src/transformations/PJ_helmert.cpp delete mode 100644 src/transformations/PJ_hgridshift.cpp delete mode 100644 src/transformations/PJ_horner.cpp delete mode 100644 src/transformations/PJ_molodensky.cpp delete mode 100644 src/transformations/PJ_vgridshift.cpp create mode 100644 src/transformations/affine.cpp create mode 100644 src/transformations/deformation.cpp create mode 100644 src/transformations/helmert.cpp create mode 100644 src/transformations/hgridshift.cpp create mode 100644 src/transformations/horner.cpp create mode 100644 src/transformations/molodensky.cpp create mode 100644 src/transformations/vgridshift.cpp (limited to 'src/transformations') diff --git a/src/transformations/PJ_affine.cpp b/src/transformations/PJ_affine.cpp deleted file mode 100644 index e2b668d3..00000000 --- a/src/transformations/PJ_affine.cpp +++ /dev/null @@ -1,250 +0,0 @@ -/************************************************************************ -* Copyright (c) 2018, Even Rouault -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* to deal in the Software without restriction, including without limitation -* the rights to use, copy, modify, merge, publish, distribute, sublicense, -* and/or sell copies of the Software, and to permit persons to whom the -* Software is furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included -* in all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -* DEALINGS IN THE SOFTWARE. -* -***********************************************************************/ -#define PJ_LIB__ - -#include -#include - -#include "proj_internal.h" -#include "proj.h" -#include "projects.h" - -PROJ_HEAD(affine, "Affine transformation"); -PROJ_HEAD(geogoffset, "Geographic Offset"); - -namespace { // anonymous namespace -struct pj_affine_coeffs { - double s11; - double s12; - double s13; - double s21; - double s22; - double s23; - double s31; - double s32; - double s33; - double tscale; -}; -} // anonymous namespace - -namespace { // anonymous namespace -struct pj_opaque_affine { - double xoff; - double yoff; - double zoff; - double toff; - struct pj_affine_coeffs forward; - struct pj_affine_coeffs reverse; -}; -} // anonymous namespace - - -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { - PJ_COORD newObs; - const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque; - const struct pj_affine_coeffs *C = &(Q->forward); - newObs.xyzt.x = Q->xoff + C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z; - newObs.xyzt.y = Q->yoff + C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z; - newObs.xyzt.z = Q->zoff + C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z; - newObs.xyzt.t = Q->toff + C->tscale * obs.xyzt.t; - return newObs; -} - -static XYZ forward_3d(LPZ lpz, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.lpz = lpz; - return forward_4d(point, P).xyz; -} - - -static XY forward_2d(LP lp, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.lp = lp; - return forward_4d(point, P).xy; -} - - -static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { - PJ_COORD newObs; - const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque; - const struct pj_affine_coeffs *C = &(Q->reverse); - obs.xyzt.x -= Q->xoff; - obs.xyzt.y -= Q->yoff; - obs.xyzt.z -= Q->zoff; - newObs.xyzt.x = C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z; - newObs.xyzt.y = C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z; - newObs.xyzt.z = C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z; - newObs.xyzt.t = C->tscale * (obs.xyzt.t - Q->toff); - return newObs; -} - -static LPZ reverse_3d(XYZ xyz, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.xyz = xyz; - return reverse_4d(point, P).lpz; -} - -static LP reverse_2d(XY xy, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.xy = xy; - return reverse_4d(point, P).lp; -} - -static struct pj_opaque_affine * initQ() { - struct pj_opaque_affine *Q = static_cast(pj_calloc(1, sizeof(struct pj_opaque_affine))); - if (nullptr==Q) - return nullptr; - - /* default values */ - Q->forward.s11 = 1.0; - Q->forward.s22 = 1.0; - Q->forward.s33 = 1.0; - Q->forward.tscale = 1.0; - - Q->reverse.s11 = 1.0; - Q->reverse.s22 = 1.0; - Q->reverse.s33 = 1.0; - Q->reverse.tscale = 1.0; - - return Q; -} - -static void computeReverseParameters(PJ* P) -{ - struct pj_opaque_affine *Q = (struct pj_opaque_affine *) P->opaque; - - /* cf https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices */ - const double a = Q->forward.s11; - const double b = Q->forward.s12; - const double c = Q->forward.s13; - const double d = Q->forward.s21; - const double e = Q->forward.s22; - const double f = Q->forward.s23; - const double g = Q->forward.s31; - const double h = Q->forward.s32; - const double i = Q->forward.s33; - const double A = e * i - f * h; - const double B = -(d * i - f * g); - const double C = (d * h - e * g); - const double D = -(b * i - c * h); - const double E = (a * i - c * g); - const double F = -(a * h - b * g); - const double G = b * f - c * e; - const double H = -(a * f - c * d); - const double I = a * e - b * d; - const double det = a * A + b * B + c * C; - if( det == 0.0 || Q->forward.tscale == 0.0 ) { - if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) { - proj_log_debug(P, "Affine: matrix non invertible"); - } - P->inv4d = nullptr; - P->inv3d = nullptr; - P->inv = nullptr; - } else { - Q->reverse.s11 = A / det; - Q->reverse.s12 = D / det; - Q->reverse.s13 = G / det; - Q->reverse.s21 = B / det; - Q->reverse.s22 = E / det; - Q->reverse.s23 = H / det; - Q->reverse.s31 = C / det; - Q->reverse.s32 = F / det; - Q->reverse.s33 = I / det; - Q->reverse.tscale = 1.0 / Q->forward.tscale; - } -} - -PJ *TRANSFORMATION(affine,0 /* no need for ellipsoid */) { - struct pj_opaque_affine *Q = initQ(); - if (nullptr==Q) - return pj_default_destructor(P, ENOMEM); - P->opaque = (void *) Q; - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = forward_2d; - P->inv = reverse_2d; - - P->left = PJ_IO_UNITS_WHATEVER; - P->right = PJ_IO_UNITS_WHATEVER; - - /* read args */ - Q->xoff = pj_param(P->ctx, P->params, "dxoff").f; - Q->yoff = pj_param(P->ctx, P->params, "dyoff").f; - Q->zoff = pj_param(P->ctx, P->params, "dzoff").f; - Q->toff = pj_param(P->ctx, P->params, "dtoff").f; - - if(pj_param (P->ctx, P->params, "ts11").i) { - Q->forward.s11 = pj_param(P->ctx, P->params, "ds11").f; - } - Q->forward.s12 = pj_param(P->ctx, P->params, "ds12").f; - Q->forward.s13 = pj_param(P->ctx, P->params, "ds13").f; - Q->forward.s21 = pj_param(P->ctx, P->params, "ds21").f; - if(pj_param (P->ctx, P->params, "ts22").i) { - Q->forward.s22 = pj_param(P->ctx, P->params, "ds22").f; - } - Q->forward.s23 = pj_param(P->ctx, P->params, "ds23").f; - Q->forward.s31 = pj_param(P->ctx, P->params, "ds31").f; - Q->forward.s32 = pj_param(P->ctx, P->params, "ds32").f; - if(pj_param (P->ctx, P->params, "ts33").i) { - Q->forward.s33 = pj_param(P->ctx, P->params, "ds33").f; - } - if(pj_param (P->ctx, P->params, "ttscale").i) { - Q->forward.tscale = pj_param(P->ctx, P->params, "dtscale").f; - } - - computeReverseParameters(P); - - return P; -} - - -/* Arcsecond to radians */ -#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) - - -PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) { - struct pj_opaque_affine *Q = initQ(); - if (nullptr==Q) - return pj_default_destructor(P, ENOMEM); - P->opaque = (void *) Q; - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = forward_2d; - P->inv = reverse_2d; - - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - - /* read args */ - Q->xoff = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD; - Q->yoff = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD; - Q->zoff = pj_param(P->ctx, P->params, "ddh").f; - - return P; -} diff --git a/src/transformations/PJ_deformation.cpp b/src/transformations/PJ_deformation.cpp deleted file mode 100644 index 6c30f21c..00000000 --- a/src/transformations/PJ_deformation.cpp +++ /dev/null @@ -1,325 +0,0 @@ -/*********************************************************************** - - Kinematic datum shifting utilizing a deformation model - - Kristian Evers, 2017-10-29 - -************************************************************************ - -Perform datum shifts by means of a deformation/velocity model. - - X_out = X_in + (T_ct - T_obs)*DX - Y_out = Y_in + (T_ct - T_obs)*DY - Z_out = Z_in + (T_ct - T_obs)*DZ - - -The deformation operation takes cartesian coordinates as input and -returns cartesian coordinates as well. - -Corrections in the gridded model are in east, north, up (ENU) space. -Hence the input coordinates needs to be converted to ENU-space when -searching for corrections in the grid. The corrections are then converted -to cartesian XYZ-space and applied to the input coordinates (also in -cartesian space). - -A full deformation model is described by two grids, one for the horizontal -components and one for the vertical component. The horizontal grid is -stored in CTable/CTable2 and the vertical grid is stored in the GTX -format. The NTv2 format should not be used for this purpose since grid- -values are scaled upon reading. Both grids are expected to contain -grid-values in units of mm/year in ENU-space. - -************************************************************************ -* Copyright (c) 2017, Kristian Evers -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* to deal in the Software without restriction, including without limitation -* the rights to use, copy, modify, merge, publish, distribute, sublicense, -* and/or sell copies of the Software, and to permit persons to whom the -* Software is furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included -* in all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -* DEALINGS IN THE SOFTWARE. -* -***********************************************************************/ -#define PJ_LIB__ -#include -#include "proj.h" -#include "proj_internal.h" -#include "proj_math.h" -#include "projects.h" - -PROJ_HEAD(deformation, "Kinematic grid shift"); - -#define TOL 1e-8 -#define MAX_ITERATIONS 10 - -namespace { // anonymous namespace -struct pj_opaque { - double t_obs; - double t_epoch; - PJ *cart; -}; -} // anonymous namespace - -/********************************************************************************/ -static XYZ get_grid_shift(PJ* P, XYZ cartesian) { -/******************************************************************************** - Read correction values from grid. The cartesian input coordinates are - converted to geodetic coordinates in order look up the correction values - in the grid. Once the grid corrections are read we need to convert them - from ENU-space to cartesian XYZ-space. ENU -> XYZ formula described in: - - Nørbech, T., et al, 2003(?), "Transformation from a Common Nordic Reference - Frame to ETRS89 in Denmark, Finland, Norway, and Sweden – status report" - -********************************************************************************/ - PJ_COORD geodetic, shift, temp; - double sp, cp, sl, cl; - int previous_errno = proj_errno_reset(P); - - /* cartesian to geodetic */ - geodetic.lpz = pj_inv3d(cartesian, static_cast(P->opaque)->cart); - - /* look up correction values in grids */ - shift.lp = proj_hgrid_value(P, geodetic.lp); - shift.enu.u = proj_vgrid_value(P, geodetic.lp); - - if (proj_errno(P) == PJD_ERR_GRID_AREA) - proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model", - proj_todeg(geodetic.lp.lam), proj_todeg(geodetic.lp.phi)); - - /* grid values are stored as mm/yr, we need m/yr */ - shift.xyz.x /= 1000; - shift.xyz.y /= 1000; - shift.xyz.z /= 1000; - - /* pre-calc cosines and sines */ - sp = sin(geodetic.lp.phi); - cp = cos(geodetic.lp.phi); - sl = sin(geodetic.lp.lam); - cl = cos(geodetic.lp.lam); - - /* ENU -> XYZ */ - temp.xyz.x = -sp*cl*shift.enu.n - sl*shift.enu.e + cp*cl*shift.enu.u; - temp.xyz.y = -sp*sl*shift.enu.n + cl*shift.enu.e + cp*sl*shift.enu.u; - temp.xyz.z = cp*shift.enu.n + sp*shift.enu.u; - - shift.xyz = temp.xyz; - - proj_errno_restore(P, previous_errno); - - return shift.xyz; -} - -/********************************************************************************/ -static XYZ reverse_shift(PJ *P, XYZ input, double dt) { -/******************************************************************************** - Iteratively determine the reverse grid shift correction values. -*********************************************************************************/ - XYZ out, delta, dif; - double z0; - int i = MAX_ITERATIONS; - - delta = get_grid_shift(P, input); - - /* Store the origial z shift for later application */ - z0 = delta.z; - - /* When iterating to find the best horizontal coordinate we also carry */ - /* along the z-component, since we need it for the cartesian -> geodetic */ - /* conversion. The z-component adjustment is overwritten with z0 after */ - /* the loop has finished. */ - out.x = input.x - dt*delta.x; - out.y = input.y - dt*delta.y; - out.z = input.z + dt*delta.z; - - do { - delta = get_grid_shift(P, out); - - if (delta.x == HUGE_VAL) - break; - - dif.x = out.x + dt*delta.x - input.x; - dif.y = out.y + dt*delta.y - input.y; - dif.z = out.z - dt*delta.z - input.z; - out.x += dif.x; - out.y += dif.y; - out.z += dif.z; - - } while ( --i && hypot(dif.x, dif.y) > TOL ); - - out.z = input.z - dt*z0; - - return out; -} - -static XYZ forward_3d(LPZ lpz, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - PJ_COORD out, in; - XYZ shift; - double dt = 0.0; - in.lpz = lpz; - out = in; - - if (Q->t_obs != HUGE_VAL) { - dt = Q->t_epoch - Q->t_obs; - } else { - out = proj_coord_error(); /* in the 3D case +t_obs must be specified */ - proj_log_debug(P, "deformation: +t_obs must be specified"); - return out.xyz; - } - - shift = get_grid_shift(P, in.xyz); - - out.xyz.x += dt * shift.x; - out.xyz.y += dt * shift.y; - out.xyz.z += dt * shift.z; - - return out.xyz; -} - - -static PJ_COORD forward_4d(PJ_COORD in, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - double dt; - XYZ shift; - PJ_COORD out = in; - - if (Q->t_obs != HUGE_VAL) { - dt = Q->t_epoch - Q->t_obs; - } else { - dt = Q->t_epoch - in.xyzt.t; - } - - shift = get_grid_shift(P, in.xyz); - - out.xyzt.x += dt*shift.x; - out.xyzt.y += dt*shift.y; - out.xyzt.z += dt*shift.z; - - - return out; -} - - -static LPZ reverse_3d(XYZ in, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - PJ_COORD out; - double dt = 0.0; - out.xyz = in; - - if (Q->t_obs != HUGE_VAL) { - dt = Q->t_epoch - Q->t_obs; - } else { - out = proj_coord_error(); /* in the 3D case +t_obs must be specified */ - proj_log_debug(P, "deformation: +t_obs must be specified"); - return out.lpz; - } - - out.xyz = reverse_shift(P, in, dt); - - return out.lpz; -} - -static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - PJ_COORD out = in; - double dt; - - - if (Q->t_obs != HUGE_VAL) { - dt = Q->t_epoch - Q->t_obs; - } else { - dt = Q->t_epoch - in.xyzt.t; - } - - out.xyz = reverse_shift(P, in.xyz, dt); - return out; -} - -static PJ *destructor(PJ *P, int errlev) { - if (nullptr==P) - return nullptr; - - if (nullptr==P->opaque) - return pj_default_destructor (P, errlev); - - if (static_cast(P->opaque)->cart) - static_cast(P->opaque)->cart->destructor (static_cast(P->opaque)->cart, errlev); - - return pj_default_destructor(P, errlev); -} - - -PJ *TRANSFORMATION(deformation,1) { - int has_xy_grids = 0; - int has_z_grids = 0; - struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); - if (nullptr==Q) - return destructor(P, ENOMEM); - P->opaque = (void *) Q; - - Q->cart = proj_create(P->ctx, "+proj=cart"); - if (Q->cart == nullptr) - return destructor(P, ENOMEM); - - /* inherit ellipsoid definition from P to Q->cart */ - pj_inherit_ellipsoid_def (P, Q->cart); - - has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i; - has_z_grids = pj_param(P->ctx, P->params, "tz_grids").i; - - /* Build gridlists. Both horizontal and vertical grids are mandatory. */ - if (!has_xy_grids || !has_z_grids) { - proj_log_error(P, "deformation: Both +xy_grids and +z_grids should be specified."); - return destructor(P, PJD_ERR_NO_ARGS ); - } - - proj_hgrid_init(P, "xy_grids"); - if (proj_errno(P)) { - proj_log_error(P, "deformation: could not find requested xy_grid(s)."); - return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); - } - - proj_vgrid_init(P, "z_grids"); - if (proj_errno(P)) { - proj_log_error(P, "deformation: could not find requested z_grid(s)."); - return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); - } - - Q->t_obs = HUGE_VAL; - if (pj_param(P->ctx, P->params, "tt_obs").i) { - Q->t_obs = pj_param(P->ctx, P->params, "dt_obs").f; - } - - if (pj_param(P->ctx, P->params, "tt_epoch").i) { - Q->t_epoch = pj_param(P->ctx, P->params, "dt_epoch").f; - } else { - proj_log_error(P, "deformation: +t_epoch parameter missing."); - return destructor(P, PJD_ERR_MISSING_ARGS); - } - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = nullptr; - P->inv = nullptr; - - P->left = PJ_IO_UNITS_CARTESIAN; - P->right = PJ_IO_UNITS_CARTESIAN; - P->destructor = destructor; - - return P; -} - diff --git a/src/transformations/PJ_helmert.cpp b/src/transformations/PJ_helmert.cpp deleted file mode 100644 index 4a3abf4e..00000000 --- a/src/transformations/PJ_helmert.cpp +++ /dev/null @@ -1,755 +0,0 @@ -/*********************************************************************** - - 3-, 4-and 7-parameter shifts, and their 6-, 8- - and 14-parameter kinematic counterparts. - - Thomas Knudsen, 2016-05-24 - -************************************************************************ - - 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. - -************************************************************************ - -Thomas Knudsen, thokn@sdfe.dk, 2016-05-24/06-05 -Kristian Evers, kreve@sdfe.dk, 2017-05-01 -Even Rouault, even.roault@spatialys.com -Last update: 2018-10-26 - -************************************************************************ -* Copyright (c) 2016, Thomas Knudsen / SDFE -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* to deal in the Software without restriction, including without limitation -* the rights to use, copy, modify, merge, publish, distribute, sublicense, -* and/or sell copies of the Software, and to permit persons to whom the -* Software is furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included -* in all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -* DEALINGS IN THE SOFTWARE. -* -***********************************************************************/ - -#define PJ_LIB__ - -#include -#include - -#include "proj_internal.h" -#include "projects.h" -#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); - - - -/***********************************************************************/ -namespace { // anonymous namespace -struct pj_opaque_helmert { -/************************************************************************ - Projection specific elements for the "helmert" PJ object -************************************************************************/ - XYZ xyz; - XYZ xyz_0; - XYZ dxyz; - XYZ refp; - PJ_OPK opk; - PJ_OPK opk_0; - PJ_OPK dopk; - double scale; - double scale_0; - double dscale; - double theta; - double theta_0; - double dtheta; - double R[3][3]; - double t_epoch, t_obs; - int no_rotation, exact, fourparam; - int is_position_vector; /* 1 = position_vector, 0 = coordinate_frame */ -}; -} // anonymous namespace - - -/* Make the maths of the rotation operations somewhat more readable and textbook like */ -#define R00 (Q->R[0][0]) -#define R01 (Q->R[0][1]) -#define R02 (Q->R[0][2]) - -#define R10 (Q->R[1][0]) -#define R11 (Q->R[1][1]) -#define R12 (Q->R[1][2]) - -#define R20 (Q->R[2][0]) -#define R21 (Q->R[2][1]) -#define R22 (Q->R[2][2]) - -/**************************************************************************/ -static void update_parameters(PJ *P) { -/*************************************************************************** - - Update transformation parameters. - --------------------------------- - - The 14-parameter Helmert transformation is at it's core the same as the - 7-parameter transformation, since the transformation parameters are - projected forward or backwards in time via the rate of changes of the - parameters. The transformation parameters are calculated for a specific - epoch before the actual Helmert transformation is carried out. - - The transformation parameters are updated with the following equation [0]: - - . - P(t) = P(EPOCH) + P * (t - EPOCH) - - . - where EPOCH is the epoch indicated in the above table and P is the rate - of that parameter. - - [0] http://itrf.ign.fr/doc_ITRF/Transfo-ITRF2008_ITRFs.txt - -*******************************************************************************/ - - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - double dt = Q->t_obs - Q->t_epoch; - - Q->xyz.x = Q->xyz_0.x + Q->dxyz.x * dt; - Q->xyz.y = Q->xyz_0.y + Q->dxyz.y * dt; - Q->xyz.z = Q->xyz_0.z + Q->dxyz.z * dt; - - Q->opk.o = Q->opk_0.o + Q->dopk.o * dt; - Q->opk.p = Q->opk_0.p + Q->dopk.p * dt; - Q->opk.k = Q->opk_0.k + Q->dopk.k * dt; - - Q->scale = Q->scale_0 + Q->dscale * dt; - - Q->theta = Q->theta_0 + Q->dtheta * dt; - - /* debugging output */ - if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) { - proj_log_trace(P, "Transformation parameters for observation " - "t_obs=%g (t_epoch=%g):", Q->t_obs, Q->t_epoch); - proj_log_trace(P, "x: %g", Q->xyz.x); - proj_log_trace(P, "y: %g", Q->xyz.y); - proj_log_trace(P, "z: %g", Q->xyz.z); - proj_log_trace(P, "s: %g", Q->scale*1e-6); - proj_log_trace(P, "rx: %g", Q->opk.o); - proj_log_trace(P, "ry: %g", Q->opk.p); - proj_log_trace(P, "rz: %g", Q->opk.k); - proj_log_trace(P, "theta: %g", Q->theta); - } -} - -/**************************************************************************/ -static void build_rot_matrix(PJ *P) { -/*************************************************************************** - - Build rotation matrix. - ---------------------- - - Here we rename rotation indices from omega, phi, kappa (opk), to - fi (i.e. phi), theta, psi (ftp), in order to reduce the mental agility - needed to implement the expression for the rotation matrix derived over - at https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions - The relevant section is Euler angles ( z-’-x" intrinsic) -> Rotation matrix - - By default small angle approximations are used: - The matrix elements are approximated by expanding the trigonometric - functions to linear order (i.e. cos(x) = 1, sin(x) = x), and discarding - products of second order. - - This was a useful hack when calculating by hand was the only option, - but in general, today, should be avoided because: - - 1. It does not save much computation time, as the rotation matrix - is built only once and probably used many times (except when - transforming spatio-temporal coordinates). - - 2. The error induced may be too large for ultra high accuracy - applications: the Earth is huge and the linear error is - approximately the angular error multiplied by the Earth radius. - - However, in many cases the approximation is necessary, since it has - been used historically: Rotation angles from older published datum - shifts may actually be a least squares fit to the linearized rotation - approximation, hence not being strictly valid for deriving the exact - rotation matrix. In fact, most publicly available transformation - parameters are based on the approximate Helmert transform, which is why - we use that as the default setting, even though it is more correct to - use the exact form of the equations. - - So in order to fit historically derived coordinates, the access to - the approximate rotation matrix is necessary - at least in principle. - - Also, when using any published datum transformation information, one - should always check which convention (exact or approximate rotation - matrix) is expected, and whether the induced error for selecting - the opposite convention is acceptable (which it often is). - - - Sign conventions - ---------------- - - Take care: Two different sign conventions exist for the rotation terms. - - Conceptually they relate to whether we rotate the coordinate system - or the "position vector" (the vector going from the coordinate system - origin to the point being transformed, i.e. the point coordinates - interpreted as vector coordinates). - - Switching between the "position vector" and "coordinate system" - conventions is simply a matter of switching the sign of the rotation - angles, which algebraically also translates into a transposition of - the rotation matrix. - - Hence, as geodetic constants should preferably be referred to exactly - as published, the "convention" option provides the ability to switch - between the conventions. - -***************************************************************************/ - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - - double f, t, p; /* phi/fi , theta, psi */ - double cf, ct, cp; /* cos (fi, theta, psi) */ - double sf, st, sp; /* sin (fi, theta, psi) */ - - /* rename (omega, phi, kappa) to (fi, theta, psi) */ - f = Q->opk.o; - t = Q->opk.p; - p = Q->opk.k; - - /* Those equations are given assuming coordinate frame convention. */ - /* For the position vector convention, we transpose the matrix just after. */ - if (Q->exact) { - cf = cos(f); - sf = sin(f); - ct = cos(t); - st = sin(t); - cp = cos(p); - sp = sin(p); - - - R00 = ct*cp; - R01 = cf*sp + sf*st*cp; - R02 = sf*sp - cf*st*cp; - - R10 = -ct*sp; - R11 = cf*cp - sf*st*sp; - R12 = sf*cp + cf*st*sp; - - R20 = st; - R21 = -sf*ct; - R22 = cf*ct; - } else{ - R00 = 1; - R01 = p; - R02 = -t; - - R10 = -p; - R11 = 1; - R12 = f; - - R20 = t; - R21 = -f; - R22 = 1; - } - - - /* - For comparison: Description from Engsager/Poder implementation - in set_dtm_1.c (trlib) - - DATUM SHIFT: - TO = scale * ROTZ * ROTY * ROTX * FROM + TRANSLA - - ( cz sz 0) (cy 0 -sy) (1 0 0) - ROTZ=(-sz cz 0), ROTY=(0 1 0), ROTX=(0 cx sx) - ( 0 0 1) (sy 0 cy) (0 -sx cx) - - trp->r11 = cos_ry*cos_rz; - trp->r12 = cos_rx*sin_rz + sin_rx*sin_ry*cos_rz; - trp->r13 = sin_rx*sin_rz - cos_rx*sin_ry*cos_rz; - - trp->r21 = -cos_ry*sin_rz; - trp->r22 = cos_rx*cos_rz - sin_rx*sin_ry*sin_rz; - trp->r23 = sin_rx*cos_rz + cos_rx*sin_ry*sin_rz; - - trp->r31 = sin_ry; - trp->r32 = -sin_rx*cos_ry; - trp->r33 = cos_rx*cos_ry; - - trp->scale = 1.0 + scale; - */ - - - if (Q->is_position_vector) { - double r; - r = R01; R01 = R10; R10 = r; - r = R02; R02 = R20; R20 = r; - r = R12; R12 = R21; R21 = r; - } - - /* some debugging output */ - if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) { - proj_log_trace(P, "Rotation Matrix:"); - proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R00, R01, R02); - proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R10, R11, R12); - proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R20, R21, R22); - } -} - - - - -/***********************************************************************/ -static XY helmert_forward (LP lp, PJ *P) { -/***********************************************************************/ - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - double x, y, cr, sr; - point.lp = lp; - - cr = cos(Q->theta) * Q->scale; - sr = sin(Q->theta) * Q->scale; - x = point.xy.x; - y = point.xy.y; - - point.xy.x = cr*x + sr*y + Q->xyz_0.x; - point.xy.y = -sr*x + cr*y + Q->xyz_0.y; - - return point.xy; -} - - -/***********************************************************************/ -static LP helmert_reverse (XY xy, PJ *P) { -/***********************************************************************/ - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - double x, y, sr, cr; - point.xy = xy; - - cr = cos(Q->theta) / Q->scale; - sr = sin(Q->theta) / Q->scale; - x = point.xy.x - Q->xyz_0.x; - y = point.xy.y - Q->xyz_0.y; - - point.xy.x = x*cr - y*sr; - point.xy.y = x*sr + y*cr; - - return point.lp; -} - - -/***********************************************************************/ -static XYZ helmert_forward_3d (LPZ lpz, PJ *P) { -/***********************************************************************/ - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - double X, Y, Z, scale; - - point.lpz = lpz; - - if (Q->fourparam) { - point.xy = helmert_forward(point.lp, P); - return point.xyz; - } - - if (Q->no_rotation) { - point.xyz.x = lpz.lam + Q->xyz.x; - point.xyz.y = lpz.phi + Q->xyz.y; - point.xyz.z = lpz.z + Q->xyz.z; - return point.xyz; - } - - scale = 1 + Q->scale * 1e-6; - - 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; /* for Molodensky-Badekas, Q->xyz already incorporates the Q->refp offset */ - point.xyz.y += Q->xyz.y; - point.xyz.z += Q->xyz.z; - - return point.xyz; -} - - -/***********************************************************************/ -static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) { -/***********************************************************************/ - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - double X, Y, Z, scale; - - point.xyz = xyz; - - if (Q->fourparam) { - point.lp = helmert_reverse(point.xy, P); - return point.lpz; - } - - if (Q->no_rotation) { - point.xyz.x = xyz.x - Q->xyz.x; - point.xyz.y = xyz.y - Q->xyz.y; - point.xyz.z = xyz.z - Q->xyz.z; - return point.lpz; - } - - scale = 1 + Q->scale * 1e-6; - - /* Unscale and deoffset */ - X = (xyz.x - Q->xyz.x) / scale; - Y = (xyz.y - Q->xyz.y) / scale; - Z = (xyz.z - Q->xyz.z) / scale; - - /* Inverse rotation through transpose multiplication */ - 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; -} - - -static PJ_COORD helmert_forward_4d (PJ_COORD point, PJ *P) { - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - - /* We only need to rebuild the rotation matrix if the - * observation time is different from the last call */ - double t_obs = (point.xyzt.t == HUGE_VAL) ? Q->t_epoch : point.xyzt.t; - if (t_obs != Q->t_obs) { - Q->t_obs = t_obs; - update_parameters(P); - build_rot_matrix(P); - } - - point.xyz = helmert_forward_3d (point.lpz, P); - - return point; -} - - -static PJ_COORD helmert_reverse_4d (PJ_COORD point, PJ *P) { - struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; - - /* We only need to rebuild the rotation matrix if the - * observation time is different from the last call */ - double t_obs = (point.xyzt.t == HUGE_VAL) ? Q->t_epoch : point.xyzt.t; - if (t_obs != Q->t_obs) { - Q->t_obs = t_obs; - update_parameters(P); - build_rot_matrix(P); - } - - point.lpz = helmert_reverse_3d (point.xyz, P); - - return point; -} - -/* Arcsecond to radians */ -#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) - - -static PJ* init_helmert_six_parameters(PJ* P) { - struct pj_opaque_helmert *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_helmert))); - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = (void *) Q; - - /* In most cases, we work on 3D cartesian coordinates */ - P->left = PJ_IO_UNITS_CARTESIAN; - P->right = PJ_IO_UNITS_CARTESIAN; - - /* 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 nullptr; - } - - /* 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. " - "Use convention=position_vector/coordinate_frame"); - return pj_default_destructor (P, PJD_ERR_INVALID_ARG); - } - - /* Support the classic PROJ towgs84 parameter, but allow later overrides.*/ - /* Note that if towgs84 is specified, the datum_params array is set up */ - /* for us automagically by the pj_datum_set call in pj_init_ctx */ - if (pj_param_exists (P->params, "towgs84")) { - Q->xyz_0.x = P->datum_params[0]; - Q->xyz_0.y = P->datum_params[1]; - Q->xyz_0.z = P->datum_params[2]; - - Q->opk_0.o = P->datum_params[3]; - Q->opk_0.p = P->datum_params[4]; - Q->opk_0.k = P->datum_params[5]; - - /* We must undo conversion to absolute scale from pj_datum_set */ - if (0==P->datum_params[6]) - Q->scale_0 = 0; - else - Q->scale_0 = (P->datum_params[6] - 1) * 1e6; - } - - 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; - Q->scale_0 = 1.0; /* default scale for the 4-param shift */ - } - - /* Scale */ - if (pj_param (P->ctx, P->params, "ts").i) { - Q->scale_0 = pj_param (P->ctx, P->params, "ds").f; - if (pj_param (P->ctx, P->params, "ttheta").i && Q->scale_0 == 0.0) - return pj_default_destructor (P, PJD_ERR_INVALID_SCALE); - } - - /* Translation rates */ - if (pj_param(P->ctx, P->params, "tdx").i) - Q->dxyz.x = pj_param (P->ctx, P->params, "ddx").f; - - if (pj_param(P->ctx, P->params, "tdy").i) - Q->dxyz.y = pj_param (P->ctx, P->params, "ddy").f; - - if (pj_param(P->ctx, P->params, "tdz").i) - Q->dxyz.z = pj_param (P->ctx, P->params, "ddz").f; - - /* Rotations rates */ - if (pj_param (P->ctx, P->params, "tdrx").i) - Q->dopk.o = pj_param (P->ctx, P->params, "ddrx").f * ARCSEC_TO_RAD; - - if (pj_param (P->ctx, P->params, "tdry").i) - Q->dopk.p = pj_param (P->ctx, P->params, "ddry").f * ARCSEC_TO_RAD; - - if (pj_param (P->ctx, P->params, "tdrz").i) - Q->dopk.k = pj_param (P->ctx, P->params, "ddrz").f * ARCSEC_TO_RAD; - - if (pj_param (P->ctx, P->params, "tdtheta").i) - Q->dtheta = pj_param (P->ctx, P->params, "ddtheta").f * ARCSEC_TO_RAD; - - /* Scale rate */ - if (pj_param (P->ctx, P->params, "tds").i) - Q->dscale = pj_param (P->ctx, P->params, "dds").f; - - - /* Epoch */ - if (pj_param(P->ctx, P->params, "tt_epoch").i) - Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f; - - if (pj_param(P->ctx, P->params, "tt_obs").i) - Q->t_obs = pj_param (P->ctx, P->params, "dt_obs").f; - - Q->xyz = Q->xyz_0; - Q->opk = Q->opk_0; - Q->scale = Q->scale_0; - Q->theta = Q->theta_0; - - if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0) && - (Q->dopk.o==0) && (Q->dopk.p==0) && (Q->dopk.k==0)) { - Q->no_rotation = 1; - } - - if( !read_convention(P) ) { - return nullptr; - } - - /* Let's help with debugging */ - if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) { - proj_log_debug(P, "Helmert parameters:"); - proj_log_debug(P, "x= %8.5f y= %8.5f z= %8.5f", Q->xyz.x, Q->xyz.y, Q->xyz.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->no_rotation ? "" : - Q->is_position_vector ? " convention=position_vector" : - " convention=coordinate_frame"); - proj_log_debug(P, "dx= %8.5f dy= %8.5f dz= %8.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z); - proj_log_debug(P, "drx=%8.5f dry=%8.5f drz=%8.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k); - proj_log_debug(P, "ds= %8.5f t_epoch=%8.5f t_obs=%8.5f", Q->dscale, Q->t_epoch, Q->t_obs); - } - - if (Q->no_rotation) { - return P; - } - - update_parameters(P); - build_rot_matrix(P); - - return P; -} - - -/***********************************************************************/ -PJ *TRANSFORMATION(molobadekas, 0) { -/***********************************************************************/ - - struct pj_opaque_helmert *Q; - - if( !init_helmert_six_parameters(P) ) { - return nullptr; - } - - 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 nullptr; - } - - /* 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/transformations/PJ_hgridshift.cpp b/src/transformations/PJ_hgridshift.cpp deleted file mode 100644 index f0e57251..00000000 --- a/src/transformations/PJ_hgridshift.cpp +++ /dev/null @@ -1,133 +0,0 @@ -#define PJ_LIB__ - -#include -#include -#include -#include - -#include "proj_internal.h" -#include "projects.h" - -PROJ_HEAD(hgridshift, "Horizontal grid shift"); - -namespace { // anonymous namespace -struct pj_opaque_hgridshift { - double t_final; - double t_epoch; -}; -} // anonymous namespace - -static XYZ forward_3d(LPZ lpz, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.lpz = lpz; - - if (P->gridlist != nullptr) { - /* Only try the gridshift if at least one grid is loaded, - * otherwise just pass the coordinate through unchanged. */ - point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD); - } - - return point.xyz; -} - - -static LPZ reverse_3d(XYZ xyz, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.xyz = xyz; - - if (P->gridlist != nullptr) { - /* Only try the gridshift if at least one grid is loaded, - * otherwise just pass the coordinate through unchanged. */ - point.lp = proj_hgrid_apply(P, point.lp, PJ_INV); - } - - return point.lpz; -} - -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque; - PJ_COORD point = obs; - - /* If transformation is not time restricted, we always call it */ - if (Q->t_final==0 || Q->t_epoch==0) { - point.xyz = forward_3d (obs.lpz, P); - return point; - } - - /* Time restricted - only apply transform if within time bracket */ - if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - point.xyz = forward_3d (obs.lpz, P); - - - return point; -} - -static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque; - PJ_COORD point = obs; - - /* If transformation is not time restricted, we always call it */ - if (Q->t_final==0 || Q->t_epoch==0) { - point.lpz = reverse_3d (obs.xyz, P); - return point; - } - - /* Time restricted - only apply transform if within time bracket */ - if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - point.lpz = reverse_3d (obs.xyz, P); - - return point; -} - - -PJ *TRANSFORMATION(hgridshift,0) { - struct pj_opaque_hgridshift *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_hgridshift))); - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = (void *) Q; - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = nullptr; - P->inv = nullptr; - - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - - if (0==pj_param(P->ctx, P->params, "tgrids").i) { - proj_log_error(P, "hgridshift: +grids parameter missing."); - return pj_default_destructor (P, PJD_ERR_NO_ARGS); - } - - /* TODO: Refactor into shared function that can be used */ - /* by both vgridshift and hgridshift */ - if (pj_param(P->ctx, P->params, "tt_final").i) { - Q->t_final = pj_param (P->ctx, P->params, "dt_final").f; - if (Q->t_final == 0) { - /* a number wasn't passed to +t_final, let's see if it was "now" */ - /* and set the time accordingly. */ - if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) { - time_t now; - struct tm *date; - time(&now); - date = localtime(&now); - Q->t_final = 1900.0 + date->tm_year + date->tm_yday/365.0; - } - } - } - - if (pj_param(P->ctx, P->params, "tt_epoch").i) - Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f; - - - proj_hgrid_init(P, "grids"); - /* Was gridlist compiled properly? */ - if ( proj_errno(P) ) { - proj_log_error(P, "hgridshift: could not find required grid(s)."); - return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); - } - - return P; -} diff --git a/src/transformations/PJ_horner.cpp b/src/transformations/PJ_horner.cpp deleted file mode 100644 index 73977de6..00000000 --- a/src/transformations/PJ_horner.cpp +++ /dev/null @@ -1,513 +0,0 @@ -/*********************************************************************** - - Interfacing to a classic piece of geodetic software - -************************************************************************ - - gen_pol is a highly efficient, classic implementation of a generic - 2D Horner's Scheme polynomial evaluation routine by Knud Poder and - Karsten Engsager, originating in the vivid geodetic environment at - what was then (1960-ish) the Danish Geodetic Institute. - - The original Poder/Engsager gen_pol implementation (where - the polynomial degree and two sets of polynomial coefficients - are packed together in one compound array, handled via a plain - double pointer) is compelling and "true to the code history": - - It has a beautiful classical 1960s ring to it, not unlike the - original fft implementations, which revolutionized spectral - analysis in twenty lines of code. - - The Poder coding sound, as classic 1960s as Phil Spector's Wall - of Sound, is beautiful and inimitable. - - On the other hand: For the uninitiated, the gen_pol code is hard - to follow, despite being compact. - - Also, since adding metadata and improving maintainability - of the code are among the implied goals of a current SDFE/DTU Space - project, the material in this file introduces a version with a - more modern (or at least 1990s) look, introducing a "double 2D - polynomial" data type, HORNER. - - Despite introducing a new data type for handling the polynomial - coefficients, great care has been taken to keep the coefficient - array organization identical to that of gen_pol. - - Hence, on one hand, the HORNER data type helps improving the - long term maintainability of the code by making the data - organization more mentally accessible. - - On the other hand, it allows us to preserve the business end of - the original gen_pol implementation - although not including the - famous "Poder dual autocheck" in all its enigmatic elegance. - - ********************************************************************** - - The material included here was written by Knud Poder, starting - around 1960, and Karsten Engsager, starting around 1970. It was - originally written in Algol 60, later (1980s) reimplemented in C. - - The HORNER data type interface, and the organization as a header - library was implemented by Thomas Knudsen, starting around 2015. - - *********************************************************************** - * - * Copyright (c) 2016, SDFE http://www.sdfe.dk / Thomas Knudsen / Karsten Engsager - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL - * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER - * DEALINGS IN THE SOFTWARE. - * - *****************************************************************************/ - -#define PJ_LIB__ - -#include -#include -#include -#include -#include - -#include "proj_internal.h" -#include "projects.h" - -PROJ_HEAD(horner, "Horner polynomial evaluation"); - -/* make horner.h interface with proj's memory management */ -#define horner_dealloc(x) pj_dealloc(x) -#define horner_calloc(n,x) pj_calloc(n,x) - -namespace { // anonymous namespace -struct horner { - int uneg; /* u axis negated? */ - int vneg; /* v axis negated? */ - int order; /* maximum degree of polynomium */ - int coefs; /* number of coefficients for each polynomium */ - double range; /* radius of the region of validity */ - - double *fwd_u; /* coefficients for the forward transformations */ - double *fwd_v; /* i.e. latitude/longitude to northing/easting */ - - double *inv_u; /* coefficients for the inverse transformations */ - double *inv_v; /* i.e. northing/easting to latitude/longitude */ - - double *fwd_c; /* coefficients for the complex forward transformations */ - double *inv_c; /* coefficients for the complex inverse transformations */ - - UV *fwd_origin; /* False longitude/latitude */ - UV *inv_origin; /* False easting/northing */ -}; -} // anonymous namespace - -typedef struct horner HORNER; -static UV horner_func (const HORNER *transformation, PJ_DIRECTION direction, UV position); -static HORNER *horner_alloc (size_t order, int complex_polynomia); -static void horner_free (HORNER *h); - -/* e.g. degree = 2: a + bx + cy + dxx + eyy + fxy, i.e. 6 coefficients */ -#define horner_number_of_coefficients(order) \ - (((order + 1)*(order + 2)/2)) - - -static void horner_free (HORNER *h) { - horner_dealloc (h->inv_v); - horner_dealloc (h->inv_u); - horner_dealloc (h->fwd_v); - horner_dealloc (h->fwd_u); - horner_dealloc (h->fwd_c); - horner_dealloc (h->inv_c); - horner_dealloc (h->fwd_origin); - horner_dealloc (h->inv_origin); - horner_dealloc (h); -} - - -static HORNER *horner_alloc (size_t order, int complex_polynomia) { - /* size_t is unsigned, so we need not check for order > 0 */ - int n = (int)horner_number_of_coefficients(order); - int polynomia_ok = 0; - HORNER *h = static_cast(horner_calloc (1, sizeof (HORNER))); - - if (nullptr==h) - return nullptr; - - if (complex_polynomia) - n = 2*(int)order + 2; - h->order = (int)order; - h->coefs = n; - - if (complex_polynomia) { - h->fwd_c = static_cast(horner_calloc (n, sizeof(double))); - h->inv_c = static_cast(horner_calloc (n, sizeof(double))); - if (h->fwd_c && h->inv_c) - polynomia_ok = 1; - } - else { - h->fwd_u = static_cast(horner_calloc (n, sizeof(double))); - h->fwd_v = static_cast(horner_calloc (n, sizeof(double))); - h->inv_u = static_cast(horner_calloc (n, sizeof(double))); - h->inv_v = static_cast(horner_calloc (n, sizeof(double))); - if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v) - polynomia_ok = 1; - } - - h->fwd_origin = static_cast(horner_calloc (1, sizeof(UV))); - h->inv_origin = static_cast(horner_calloc (1, sizeof(UV))); - - if (polynomia_ok && h->fwd_origin && h->inv_origin) - return h; - - /* safe, since all pointers are null-initialized (by calloc) */ - horner_free (h); - return nullptr; -} - - - - -/**********************************************************************/ -static UV horner_func (const HORNER *transformation, PJ_DIRECTION direction, UV position) { -/*********************************************************************** - -A reimplementation of the classic Engsager/Poder 2D Horner polynomial -evaluation engine "gen_pol". - -This version omits the inimitable Poder "dual autocheck"-machinery, -which here is intended to be implemented at a higher level of the -library: We separate the polynomial evaluation from the quality -control (which, given the limited MTBF for "computing machinery", -typical when Knud Poder invented the dual autocheck method, -was not defensible at that time). - -Another difference from the original version is that we return the -result on the stack, rather than accepting pointers to result variables -as input. This results in code that is easy to read: - - projected = horner (s34j, 1, geographic); - geographic = horner (s34j, -1, projected ); - -and experiments have shown that on contemporary architectures, the time -taken for returning even comparatively large objects on the stack (and -the UV is not that large - typically only 16 bytes) is negligibly -different from passing two pointers (i.e. typically also 16 bytes) the -other way. - -The polynomium has the form: - -P = sum (i = [0 : order]) - sum (j = [0 : order - i]) - pow(par_1, i) * pow(par_2, j) * coef(index(order, i, j)) - -For numerical stability, the summation is carried out backwards, -summing the tiny high order elements first. - -***********************************************************************/ - - /* These variable names follow the Engsager/Poder implementation */ - int sz; /* Number of coefficients per polynomial */ - double *tcx, *tcy; /* Coefficient pointers */ - double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ - double n, e; - UV uv_error; - uv_error.u = uv_error.v = HUGE_VAL; - - if (nullptr==transformation) - return uv_error; - - /* Check for valid value of direction (-1, 0, 1) */ - switch (direction) { - case PJ_IDENT: /* no-op */ - return position; - case PJ_FWD: /* forward */ - case PJ_INV: /* inverse */ - break; - default: /* invalid */ - errno = EINVAL; - return uv_error; - } - - /* Prepare for double Horner */ - sz = horner_number_of_coefficients(transformation->order); - range = transformation->range; - - - if (direction==PJ_FWD) { /* forward */ - tcx = transformation->fwd_u + sz; - tcy = transformation->fwd_v + sz; - e = position.u - transformation->fwd_origin->u; - n = position.v - transformation->fwd_origin->v; - } else { /* inverse */ - tcx = transformation->inv_u + sz; - tcy = transformation->inv_v + sz; - e = position.u - transformation->inv_origin->u; - n = position.v - transformation->inv_origin->v; - } - - if ((fabs(n) > range) || (fabs(e) > range)) { - errno = EDOM; - return uv_error; - } - - /* The melody of this block is straight out of the great Engsager/Poder songbook */ - else { - int g = transformation->order; - int r = g, c; - double u, v, N, E; - - /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */ - N = *--tcy; - E = *--tcx; - for (; r > 0; r--) { - u = *--tcy; - v = *--tcx; - for (c = g; c >= r; c--) { - u = n*u + *--tcy; - v = e*v + *--tcx; - } - N = e*N + u; - E = n*E + v; - } - - position.u = E; - position.v = N; - } - - return position; -} - - - - - - - -static PJ_COORD horner_forward_4d (PJ_COORD point, PJ *P) { - point.uv = horner_func ((HORNER *) P->opaque, PJ_FWD, point.uv); - return point; -} - -static PJ_COORD horner_reverse_4d (PJ_COORD point, PJ *P) { - point.uv = horner_func ((HORNER *) P->opaque, PJ_INV, point.uv); - return point; -} - - - - -/**********************************************************************/ -static UV complex_horner (const HORNER *transformation, PJ_DIRECTION direction, UV position) { -/*********************************************************************** - -A reimplementation of a classic Engsager/Poder Horner complex -polynomial evaluation engine. - -***********************************************************************/ - - /* These variable names follow the Engsager/Poder implementation */ - int sz; /* Number of coefficients */ - double *c, *cb; /* Coefficient pointers */ - double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ - double n, e, w, N, E; - UV uv_error; - uv_error.u = uv_error.v = HUGE_VAL; - - if (nullptr==transformation) - return uv_error; - - /* Check for valid value of direction (-1, 0, 1) */ - switch (direction) { - case PJ_IDENT: /* no-op */ - return position; - case PJ_FWD: /* forward */ - case PJ_INV: /* inverse */ - break; - default: /* invalid */ - errno = EINVAL; - return uv_error; - } - - /* Prepare for double Horner */ - sz = 2*transformation->order + 2; - range = transformation->range; - - if (direction==PJ_FWD) { /* forward */ - cb = transformation->fwd_c; - c = cb + sz; - e = position.u - transformation->fwd_origin->u; - n = position.v - transformation->fwd_origin->v; - if (transformation->uneg) - e = -e; - if (transformation->vneg) - n = -n; - } else { /* inverse */ - cb = transformation->inv_c; - c = cb + sz; - e = position.u - transformation->inv_origin->u; - n = position.v - transformation->inv_origin->v; - if (transformation->uneg) - e = -e; - if (transformation->vneg) - n = -n; - } - - if ((fabs(n) > range) || (fabs(e) > range)) { - errno = EDOM; - return uv_error; - } - - /* Everything's set up properly - now do the actual polynomium evaluation */ - E = *--c; - N = *--c; - while (c > cb) { - w = n*E + e*N + *--c; - N = n*N - e*E + *--c; - E = w; - } - - position.u = E; - position.v = N; - return position; -} - - - -static PJ_COORD complex_horner_forward_4d (PJ_COORD point, PJ *P) { - point.uv = complex_horner ((HORNER *) P->opaque, PJ_FWD, point.uv); - return point; -} - -static PJ_COORD complex_horner_reverse_4d (PJ_COORD point, PJ *P) { - point.uv = complex_horner ((HORNER *) P->opaque, PJ_INV, point.uv); - return point; -} - - -static PJ *horner_freeup (PJ *P, int errlev) { /* Destructor */ - if (nullptr==P) - return nullptr; - if (nullptr==P->opaque) - return pj_default_destructor (P, errlev); - horner_free ((HORNER *) P->opaque); - P->opaque = nullptr; - return pj_default_destructor (P, errlev); -} - - -static int parse_coefs (PJ *P, double *coefs, const char *param, int ncoefs) { - char *buf, *init, *next = nullptr; - int i; - - buf = static_cast(pj_calloc (strlen (param) + 2, sizeof(char))); - if (nullptr==buf) { - proj_log_error (P, "Horner: No memory left"); - return 0; - } - - sprintf (buf, "t%s", param); - if (0==pj_param (P->ctx, P->params, buf).i) { - pj_dealloc (buf); - return 0; - } - sprintf (buf, "s%s", param); - init = pj_param(P->ctx, P->params, buf).s; - pj_dealloc (buf); - - for (i = 0; i < ncoefs; i++) { - if (i > 0) { - if ( next == nullptr || ','!=*next) { - proj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs); - return 0; - } - init = ++next; - } - coefs[i] = pj_strtod (init, &next); - } - return 1; -} - - -/*********************************************************************/ -PJ *PROJECTION(horner) { -/*********************************************************************/ - int degree = 0, n, complex_polynomia = 0; - HORNER *Q; - P->fwd4d = horner_forward_4d; - P->inv4d = horner_reverse_4d; - P->fwd3d = nullptr; - P->inv3d = nullptr; - P->fwd = nullptr; - P->inv = nullptr; - P->left = P->right = PJ_IO_UNITS_PROJECTED; - P->destructor = horner_freeup; - - /* Polynomial degree specified? */ - if (pj_param (P->ctx, P->params, "tdeg").i) { /* degree specified? */ - degree = pj_param(P->ctx, P->params, "ideg").i; - if (degree < 0 || degree > 10000) { - /* What are reasonable minimum and maximums for degree? */ - proj_log_debug (P, "Horner: Degree is unreasonable: %d", degree); - return horner_freeup (P, PJD_ERR_INVALID_ARG); - } - } else { - proj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)"); - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - } - - if (pj_param (P->ctx, P->params, "tfwd_c").i || pj_param (P->ctx, P->params, "tinv_c").i) /* complex polynomium? */ - complex_polynomia = 1; - - Q = horner_alloc (degree, complex_polynomia); - if (Q == nullptr) - return horner_freeup (P, ENOMEM); - P->opaque = Q; - - if (complex_polynomia) { - /* Westings and/or southings? */ - Q->uneg = pj_param_exists (P->params, "uneg") ? 1 : 0; - Q->vneg = pj_param_exists (P->params, "vneg") ? 1 : 0; - - n = 2*degree + 2; - if (0==parse_coefs (P, Q->fwd_c, "fwd_c", n)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - if (0==parse_coefs (P, Q->inv_c, "inv_c", n)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - P->fwd4d = complex_horner_forward_4d; - P->inv4d = complex_horner_reverse_4d; - } - - else { - n = horner_number_of_coefficients (degree); - if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - if (0==parse_coefs (P, Q->inv_u, "inv_u", n)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - if (0==parse_coefs (P, Q->inv_v, "inv_v", n)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - } - - if (0==parse_coefs (P, (double *)(Q->fwd_origin), "fwd_origin", 2)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) - return horner_freeup (P, PJD_ERR_MISSING_ARGS); - if (0==parse_coefs (P, &Q->range, "range", 1)) - Q->range = 500000; - - return P; -} diff --git a/src/transformations/PJ_molodensky.cpp b/src/transformations/PJ_molodensky.cpp deleted file mode 100644 index 91743fda..00000000 --- a/src/transformations/PJ_molodensky.cpp +++ /dev/null @@ -1,329 +0,0 @@ -/*********************************************************************** - - (Abridged) Molodensky Transform - - Kristian Evers, 2017-07-07 - -************************************************************************ - - Implements the (abridged) Molodensky transformations for 2D and 3D - data. - - Primarily useful for implementation of datum shifts in transformation - pipelines. - - The code in this file is mostly based on - - The Standard and Abridged Molodensky Coordinate Transformation - Formulae, 2004, R.E. Deaking, - http://www.mygeodesy.id.au/documents/Molodensky%20V2.pdf - - - -************************************************************************ -* Copyright (c) 2017, Kristian Evers / SDFE -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* to deal in the Software without restriction, including without limitation -* the rights to use, copy, modify, merge, publish, distribute, sublicense, -* and/or sell copies of the Software, and to permit persons to whom the -* Software is furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included -* in all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -* DEALINGS IN THE SOFTWARE. -* -***********************************************************************/ -#define PJ_LIB__ - -#include -#include - -#include "proj.h" -#include "proj_internal.h" -#include "projects.h" - -PROJ_HEAD(molodensky, "Molodensky transform"); - -static XYZ forward_3d(LPZ lpz, PJ *P); -static LPZ reverse_3d(XYZ xyz, PJ *P); - -namespace { // anonymous namespace -struct pj_opaque_molodensky { - double dx; - double dy; - double dz; - double da; - double df; - int abridged; -}; -} // anonymous namespace - - -static double RN (double a, double es, double phi) { -/********************************************************** - N(phi) - prime vertical radius of curvature - ------------------------------------------- - - This is basically the same function as in PJ_cart.c - should probably be refactored into it's own file at some - point. - -**********************************************************/ - double s = sin(phi); - if (es==0) - return a; - - return a / sqrt (1 - es*s*s); -} - - -static double RM (double a, double es, double phi) { -/********************************************************** - M(phi) - Meridian radius of curvature - ------------------------------------- - - Source: - - E.J Krakiwsky & D.B. Thomson, 1974, - GEODETIC POSITION COMPUTATIONS, - - Fredericton NB, Canada: - University of New Brunswick, - Department of Geodesy and Geomatics Engineering, - Lecture Notes No. 39, - 99 pp. - - http://www2.unb.ca/gge/Pubs/LN39.pdf - -**********************************************************/ - double s = sin(phi); - if (es==0) - return a; - - /* eq. 13a */ - if (phi == 0) - return a * (1-es); - - /* eq. 13b */ - if (fabs(phi) == M_PI_2) - return a / sqrt(1-es); - - /* eq. 13 */ - return (a * (1 - es) ) / pow(1 - es*s*s, 1.5); - -} - - -static LPZ calc_standard_params(LPZ lpz, PJ *P) { - struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; - double dphi, dlam, dh; - - /* sines and cosines */ - double slam = sin(lpz.lam); - double clam = cos(lpz.lam); - double sphi = sin(lpz.phi); - double cphi = cos(lpz.phi); - - /* ellipsoid parameters and differences */ - double f = P->f, a = P->a; - double dx = Q->dx, dy = Q->dy, dz = Q->dz; - double da = Q->da, df = Q->df; - - /* ellipsoid radii of curvature */ - double rho = RM(a, P->es, lpz.phi); - double nu = RN(a, P->e2s, lpz.phi); - - /* delta phi */ - dphi = (-dx*sphi*clam) - (dy*sphi*slam) + (dz*cphi) - + ((nu * P->es * sphi * cphi * da) / a) - + (sphi*cphi * ( rho/(1-f) + nu*(1-f))*df); - dphi /= (rho + lpz.z); - - /* delta lambda */ - dlam = (-dx*slam + dy*clam) / ((nu+lpz.z)*cphi); - - /* delta h */ - dh = dx*cphi*clam + dy*cphi*slam + dz*sphi - (a/nu)*da + nu*(1-f)*sphi*sphi*df; - - lpz.phi = dphi; - lpz.lam = dlam; - lpz.z = dh; - - return lpz; -} - - -static LPZ calc_abridged_params(LPZ lpz, PJ *P) { - struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; - double dphi, dlam, dh; - - /* sines and cosines */ - double slam = sin(lpz.lam); - double clam = cos(lpz.lam); - double sphi = sin(lpz.phi); - double cphi = cos(lpz.phi); - - /* ellipsoid parameters and differences */ - double dx = Q->dx, dy = Q->dy, dz = Q->dz; - double da = Q->da, df = Q->df; - double adffda = (P->a*df + P->f*da); - - /* delta phi */ - dphi = -dx*sphi*clam - dy*sphi*slam + dz*cphi + adffda*sin(2*lpz.phi); - dphi /= RM(P->a, P->es, lpz.phi); - - /* delta lambda */ - dlam = -dx*slam + dy*clam; - dlam /= RN(P->a, P->e2s, lpz.phi)*cphi; - - /* delta h */ - dh = dx*cphi*clam + dy*cphi*slam + dz*sphi - da + adffda*sphi*sphi; - - /* offset coordinate */ - lpz.phi = dphi; - lpz.lam = dlam; - lpz.z = dh; - - return lpz; -} - - -static XY forward_2d(LP lp, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - - point.lp = lp; - point.xyz = forward_3d(point.lpz, P); - - return point.xy; -} - - -static LP reverse_2d(XY xy, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - - point.xy = xy; - point.xyz.z = 0; - point.lpz = reverse_3d(point.xyz, P); - - return point.lp; -} - - -static XYZ forward_3d(LPZ lpz, PJ *P) { - struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - - point.lpz = lpz; - - /* calculate parameters depending on the mode we are in */ - if (Q->abridged) { - lpz = calc_abridged_params(lpz, P); - } else { - lpz = calc_standard_params(lpz, P); - } - - /* offset coordinate */ - point.lpz.phi += lpz.phi; - point.lpz.lam += lpz.lam; - point.lpz.z += lpz.z; - - return point.xyz; -} - - -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { - obs.xyz = forward_3d(obs.lpz, P); - return obs; -} - - -static LPZ reverse_3d(XYZ xyz, PJ *P) { - struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - LPZ lpz; - - /* calculate parameters depending on the mode we are in */ - point.xyz = xyz; - if (Q->abridged) - lpz = calc_abridged_params(point.lpz, P); - else - lpz = calc_standard_params(point.lpz, P); - - /* offset coordinate */ - point.lpz.phi -= lpz.phi; - point.lpz.lam -= lpz.lam; - point.lpz.z -= lpz.z; - - return point.lpz; -} - - -static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { - obs.lpz = reverse_3d(obs.xyz, P); - return obs; -} - - -PJ *TRANSFORMATION(molodensky,1) { - int count_required_params = 0; - struct pj_opaque_molodensky *Q = static_cast(pj_calloc(1, sizeof(struct pj_opaque_molodensky))); - if (nullptr==Q) - return pj_default_destructor(P, ENOMEM); - P->opaque = (void *) Q; - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = forward_2d; - P->inv = reverse_2d; - - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - - /* read args */ - if (pj_param(P->ctx, P->params, "tdx").i) { - count_required_params ++; - Q->dx = pj_param(P->ctx, P->params, "ddx").f; - } - - if (pj_param(P->ctx, P->params, "tdy").i) { - count_required_params ++; - Q->dy = pj_param(P->ctx, P->params, "ddy").f; - } - - if (pj_param(P->ctx, P->params, "tdz").i) { - count_required_params ++; - Q->dz = pj_param(P->ctx, P->params, "ddz").f; - } - - if (pj_param(P->ctx, P->params, "tda").i) { - count_required_params ++; - Q->da = pj_param(P->ctx, P->params, "dda").f; - } - - if (pj_param(P->ctx, P->params, "tdf").i) { - count_required_params ++; - Q->df = pj_param(P->ctx, P->params, "ddf").f; - } - - Q->abridged = pj_param(P->ctx, P->params, "tabridged").i; - - /* We want all parameters (except +abridged) to be set */ - if (count_required_params == 0) - return pj_default_destructor(P, PJD_ERR_NO_ARGS); - - if (count_required_params != 5) - return pj_default_destructor(P, PJD_ERR_MISSING_ARGS); - - return P; -} diff --git a/src/transformations/PJ_vgridshift.cpp b/src/transformations/PJ_vgridshift.cpp deleted file mode 100644 index b3da906d..00000000 --- a/src/transformations/PJ_vgridshift.cpp +++ /dev/null @@ -1,144 +0,0 @@ -#define PJ_LIB__ - -#include -#include -#include -#include - -#include "proj_internal.h" -#include "projects.h" - -PROJ_HEAD(vgridshift, "Vertical grid shift"); - -namespace { // anonymous namespace -struct pj_opaque_vgridshift { - double t_final; - double t_epoch; - double forward_multiplier; -}; -} // anonymous namespace - -static XYZ forward_3d(LPZ lpz, PJ *P) { - struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - point.lpz = lpz; - - if (P->vgridlist_geoid != nullptr) { - /* Only try the gridshift if at least one grid is loaded, - * otherwise just pass the coordinate through unchanged. */ - point.xyz.z += Q->forward_multiplier * proj_vgrid_value(P, point.lp); - } - - return point.xyz; -} - - -static LPZ reverse_3d(XYZ xyz, PJ *P) { - struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - point.xyz = xyz; - - if (P->vgridlist_geoid != nullptr) { - /* Only try the gridshift if at least one grid is loaded, - * otherwise just pass the coordinate through unchanged. */ - point.xyz.z -= Q->forward_multiplier * proj_vgrid_value(P, point.lp); - } - - return point.lpz; -} - - -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; - PJ_COORD point = obs; - - /* If transformation is not time restricted, we always call it */ - if (Q->t_final==0 || Q->t_epoch==0) { - point.xyz = forward_3d (obs.lpz, P); - return point; - } - - /* Time restricted - only apply transform if within time bracket */ - if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - point.xyz = forward_3d (obs.lpz, P); - - - return point; -} - -static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; - PJ_COORD point = obs; - - /* If transformation is not time restricted, we always call it */ - if (Q->t_final==0 || Q->t_epoch==0) { - point.lpz = reverse_3d (obs.xyz, P); - return point; - } - - /* Time restricted - only apply transform if within time bracket */ - if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - point.lpz = reverse_3d (obs.xyz, P); - - return point; -} - - -PJ *TRANSFORMATION(vgridshift,0) { - struct pj_opaque_vgridshift *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_vgridshift))); - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = (void *) Q; - - if (!pj_param(P->ctx, P->params, "tgrids").i) { - proj_log_error(P, "vgridshift: +grids parameter missing."); - return pj_default_destructor(P, PJD_ERR_NO_ARGS); - } - - /* TODO: Refactor into shared function that can be used */ - /* by both vgridshift and hgridshift */ - if (pj_param(P->ctx, P->params, "tt_final").i) { - Q->t_final = pj_param (P->ctx, P->params, "dt_final").f; - if (Q->t_final == 0) { - /* a number wasn't passed to +t_final, let's see if it was "now" */ - /* and set the time accordingly. */ - if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) { - time_t now; - struct tm *date; - time(&now); - date = localtime(&now); - Q->t_final = 1900.0 + date->tm_year + date->tm_yday/365.0; - } - } - } - - if (pj_param(P->ctx, P->params, "tt_epoch").i) - Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f; - - /* historical: the forward direction subtracts the grid offset. */ - Q->forward_multiplier = -1.0; - if (pj_param(P->ctx, P->params, "tmultiplier").i) { - Q->forward_multiplier = pj_param(P->ctx, P->params, "dmultiplier").f; - } - - /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */ - proj_vgrid_init(P, "grids"); - - /* Was gridlist compiled properly? */ - if ( proj_errno(P) ) { - proj_log_error(P, "vgridshift: could not find required grid(s)."); - return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); - } - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = nullptr; - P->inv = nullptr; - - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - - return P; -} diff --git a/src/transformations/affine.cpp b/src/transformations/affine.cpp new file mode 100644 index 00000000..e2b668d3 --- /dev/null +++ b/src/transformations/affine.cpp @@ -0,0 +1,250 @@ +/************************************************************************ +* Copyright (c) 2018, Even Rouault +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +* DEALINGS IN THE SOFTWARE. +* +***********************************************************************/ +#define PJ_LIB__ + +#include +#include + +#include "proj_internal.h" +#include "proj.h" +#include "projects.h" + +PROJ_HEAD(affine, "Affine transformation"); +PROJ_HEAD(geogoffset, "Geographic Offset"); + +namespace { // anonymous namespace +struct pj_affine_coeffs { + double s11; + double s12; + double s13; + double s21; + double s22; + double s23; + double s31; + double s32; + double s33; + double tscale; +}; +} // anonymous namespace + +namespace { // anonymous namespace +struct pj_opaque_affine { + double xoff; + double yoff; + double zoff; + double toff; + struct pj_affine_coeffs forward; + struct pj_affine_coeffs reverse; +}; +} // anonymous namespace + + +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { + PJ_COORD newObs; + const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque; + const struct pj_affine_coeffs *C = &(Q->forward); + newObs.xyzt.x = Q->xoff + C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z; + newObs.xyzt.y = Q->yoff + C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z; + newObs.xyzt.z = Q->zoff + C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z; + newObs.xyzt.t = Q->toff + C->tscale * obs.xyzt.t; + return newObs; +} + +static XYZ forward_3d(LPZ lpz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lpz = lpz; + return forward_4d(point, P).xyz; +} + + +static XY forward_2d(LP lp, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lp = lp; + return forward_4d(point, P).xy; +} + + +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { + PJ_COORD newObs; + const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque; + const struct pj_affine_coeffs *C = &(Q->reverse); + obs.xyzt.x -= Q->xoff; + obs.xyzt.y -= Q->yoff; + obs.xyzt.z -= Q->zoff; + newObs.xyzt.x = C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z; + newObs.xyzt.y = C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z; + newObs.xyzt.z = C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z; + newObs.xyzt.t = C->tscale * (obs.xyzt.t - Q->toff); + return newObs; +} + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xyz = xyz; + return reverse_4d(point, P).lpz; +} + +static LP reverse_2d(XY xy, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xy = xy; + return reverse_4d(point, P).lp; +} + +static struct pj_opaque_affine * initQ() { + struct pj_opaque_affine *Q = static_cast(pj_calloc(1, sizeof(struct pj_opaque_affine))); + if (nullptr==Q) + return nullptr; + + /* default values */ + Q->forward.s11 = 1.0; + Q->forward.s22 = 1.0; + Q->forward.s33 = 1.0; + Q->forward.tscale = 1.0; + + Q->reverse.s11 = 1.0; + Q->reverse.s22 = 1.0; + Q->reverse.s33 = 1.0; + Q->reverse.tscale = 1.0; + + return Q; +} + +static void computeReverseParameters(PJ* P) +{ + struct pj_opaque_affine *Q = (struct pj_opaque_affine *) P->opaque; + + /* cf https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices */ + const double a = Q->forward.s11; + const double b = Q->forward.s12; + const double c = Q->forward.s13; + const double d = Q->forward.s21; + const double e = Q->forward.s22; + const double f = Q->forward.s23; + const double g = Q->forward.s31; + const double h = Q->forward.s32; + const double i = Q->forward.s33; + const double A = e * i - f * h; + const double B = -(d * i - f * g); + const double C = (d * h - e * g); + const double D = -(b * i - c * h); + const double E = (a * i - c * g); + const double F = -(a * h - b * g); + const double G = b * f - c * e; + const double H = -(a * f - c * d); + const double I = a * e - b * d; + const double det = a * A + b * B + c * C; + if( det == 0.0 || Q->forward.tscale == 0.0 ) { + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) { + proj_log_debug(P, "Affine: matrix non invertible"); + } + P->inv4d = nullptr; + P->inv3d = nullptr; + P->inv = nullptr; + } else { + Q->reverse.s11 = A / det; + Q->reverse.s12 = D / det; + Q->reverse.s13 = G / det; + Q->reverse.s21 = B / det; + Q->reverse.s22 = E / det; + Q->reverse.s23 = H / det; + Q->reverse.s31 = C / det; + Q->reverse.s32 = F / det; + Q->reverse.s33 = I / det; + Q->reverse.tscale = 1.0 / Q->forward.tscale; + } +} + +PJ *TRANSFORMATION(affine,0 /* no need for ellipsoid */) { + struct pj_opaque_affine *Q = initQ(); + if (nullptr==Q) + return pj_default_destructor(P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_WHATEVER; + P->right = PJ_IO_UNITS_WHATEVER; + + /* read args */ + Q->xoff = pj_param(P->ctx, P->params, "dxoff").f; + Q->yoff = pj_param(P->ctx, P->params, "dyoff").f; + Q->zoff = pj_param(P->ctx, P->params, "dzoff").f; + Q->toff = pj_param(P->ctx, P->params, "dtoff").f; + + if(pj_param (P->ctx, P->params, "ts11").i) { + Q->forward.s11 = pj_param(P->ctx, P->params, "ds11").f; + } + Q->forward.s12 = pj_param(P->ctx, P->params, "ds12").f; + Q->forward.s13 = pj_param(P->ctx, P->params, "ds13").f; + Q->forward.s21 = pj_param(P->ctx, P->params, "ds21").f; + if(pj_param (P->ctx, P->params, "ts22").i) { + Q->forward.s22 = pj_param(P->ctx, P->params, "ds22").f; + } + Q->forward.s23 = pj_param(P->ctx, P->params, "ds23").f; + Q->forward.s31 = pj_param(P->ctx, P->params, "ds31").f; + Q->forward.s32 = pj_param(P->ctx, P->params, "ds32").f; + if(pj_param (P->ctx, P->params, "ts33").i) { + Q->forward.s33 = pj_param(P->ctx, P->params, "ds33").f; + } + if(pj_param (P->ctx, P->params, "ttscale").i) { + Q->forward.tscale = pj_param(P->ctx, P->params, "dtscale").f; + } + + computeReverseParameters(P); + + return P; +} + + +/* Arcsecond to radians */ +#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) + + +PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) { + struct pj_opaque_affine *Q = initQ(); + if (nullptr==Q) + return pj_default_destructor(P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + /* read args */ + Q->xoff = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD; + Q->yoff = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD; + Q->zoff = pj_param(P->ctx, P->params, "ddh").f; + + return P; +} diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp new file mode 100644 index 00000000..6c30f21c --- /dev/null +++ b/src/transformations/deformation.cpp @@ -0,0 +1,325 @@ +/*********************************************************************** + + Kinematic datum shifting utilizing a deformation model + + Kristian Evers, 2017-10-29 + +************************************************************************ + +Perform datum shifts by means of a deformation/velocity model. + + X_out = X_in + (T_ct - T_obs)*DX + Y_out = Y_in + (T_ct - T_obs)*DY + Z_out = Z_in + (T_ct - T_obs)*DZ + + +The deformation operation takes cartesian coordinates as input and +returns cartesian coordinates as well. + +Corrections in the gridded model are in east, north, up (ENU) space. +Hence the input coordinates needs to be converted to ENU-space when +searching for corrections in the grid. The corrections are then converted +to cartesian XYZ-space and applied to the input coordinates (also in +cartesian space). + +A full deformation model is described by two grids, one for the horizontal +components and one for the vertical component. The horizontal grid is +stored in CTable/CTable2 and the vertical grid is stored in the GTX +format. The NTv2 format should not be used for this purpose since grid- +values are scaled upon reading. Both grids are expected to contain +grid-values in units of mm/year in ENU-space. + +************************************************************************ +* Copyright (c) 2017, Kristian Evers +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +* DEALINGS IN THE SOFTWARE. +* +***********************************************************************/ +#define PJ_LIB__ +#include +#include "proj.h" +#include "proj_internal.h" +#include "proj_math.h" +#include "projects.h" + +PROJ_HEAD(deformation, "Kinematic grid shift"); + +#define TOL 1e-8 +#define MAX_ITERATIONS 10 + +namespace { // anonymous namespace +struct pj_opaque { + double t_obs; + double t_epoch; + PJ *cart; +}; +} // anonymous namespace + +/********************************************************************************/ +static XYZ get_grid_shift(PJ* P, XYZ cartesian) { +/******************************************************************************** + Read correction values from grid. The cartesian input coordinates are + converted to geodetic coordinates in order look up the correction values + in the grid. Once the grid corrections are read we need to convert them + from ENU-space to cartesian XYZ-space. ENU -> XYZ formula described in: + + Nørbech, T., et al, 2003(?), "Transformation from a Common Nordic Reference + Frame to ETRS89 in Denmark, Finland, Norway, and Sweden – status report" + +********************************************************************************/ + PJ_COORD geodetic, shift, temp; + double sp, cp, sl, cl; + int previous_errno = proj_errno_reset(P); + + /* cartesian to geodetic */ + geodetic.lpz = pj_inv3d(cartesian, static_cast(P->opaque)->cart); + + /* look up correction values in grids */ + shift.lp = proj_hgrid_value(P, geodetic.lp); + shift.enu.u = proj_vgrid_value(P, geodetic.lp); + + if (proj_errno(P) == PJD_ERR_GRID_AREA) + proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model", + proj_todeg(geodetic.lp.lam), proj_todeg(geodetic.lp.phi)); + + /* grid values are stored as mm/yr, we need m/yr */ + shift.xyz.x /= 1000; + shift.xyz.y /= 1000; + shift.xyz.z /= 1000; + + /* pre-calc cosines and sines */ + sp = sin(geodetic.lp.phi); + cp = cos(geodetic.lp.phi); + sl = sin(geodetic.lp.lam); + cl = cos(geodetic.lp.lam); + + /* ENU -> XYZ */ + temp.xyz.x = -sp*cl*shift.enu.n - sl*shift.enu.e + cp*cl*shift.enu.u; + temp.xyz.y = -sp*sl*shift.enu.n + cl*shift.enu.e + cp*sl*shift.enu.u; + temp.xyz.z = cp*shift.enu.n + sp*shift.enu.u; + + shift.xyz = temp.xyz; + + proj_errno_restore(P, previous_errno); + + return shift.xyz; +} + +/********************************************************************************/ +static XYZ reverse_shift(PJ *P, XYZ input, double dt) { +/******************************************************************************** + Iteratively determine the reverse grid shift correction values. +*********************************************************************************/ + XYZ out, delta, dif; + double z0; + int i = MAX_ITERATIONS; + + delta = get_grid_shift(P, input); + + /* Store the origial z shift for later application */ + z0 = delta.z; + + /* When iterating to find the best horizontal coordinate we also carry */ + /* along the z-component, since we need it for the cartesian -> geodetic */ + /* conversion. The z-component adjustment is overwritten with z0 after */ + /* the loop has finished. */ + out.x = input.x - dt*delta.x; + out.y = input.y - dt*delta.y; + out.z = input.z + dt*delta.z; + + do { + delta = get_grid_shift(P, out); + + if (delta.x == HUGE_VAL) + break; + + dif.x = out.x + dt*delta.x - input.x; + dif.y = out.y + dt*delta.y - input.y; + dif.z = out.z - dt*delta.z - input.z; + out.x += dif.x; + out.y += dif.y; + out.z += dif.z; + + } while ( --i && hypot(dif.x, dif.y) > TOL ); + + out.z = input.z - dt*z0; + + return out; +} + +static XYZ forward_3d(LPZ lpz, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + PJ_COORD out, in; + XYZ shift; + double dt = 0.0; + in.lpz = lpz; + out = in; + + if (Q->t_obs != HUGE_VAL) { + dt = Q->t_epoch - Q->t_obs; + } else { + out = proj_coord_error(); /* in the 3D case +t_obs must be specified */ + proj_log_debug(P, "deformation: +t_obs must be specified"); + return out.xyz; + } + + shift = get_grid_shift(P, in.xyz); + + out.xyz.x += dt * shift.x; + out.xyz.y += dt * shift.y; + out.xyz.z += dt * shift.z; + + return out.xyz; +} + + +static PJ_COORD forward_4d(PJ_COORD in, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + double dt; + XYZ shift; + PJ_COORD out = in; + + if (Q->t_obs != HUGE_VAL) { + dt = Q->t_epoch - Q->t_obs; + } else { + dt = Q->t_epoch - in.xyzt.t; + } + + shift = get_grid_shift(P, in.xyz); + + out.xyzt.x += dt*shift.x; + out.xyzt.y += dt*shift.y; + out.xyzt.z += dt*shift.z; + + + return out; +} + + +static LPZ reverse_3d(XYZ in, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + PJ_COORD out; + double dt = 0.0; + out.xyz = in; + + if (Q->t_obs != HUGE_VAL) { + dt = Q->t_epoch - Q->t_obs; + } else { + out = proj_coord_error(); /* in the 3D case +t_obs must be specified */ + proj_log_debug(P, "deformation: +t_obs must be specified"); + return out.lpz; + } + + out.xyz = reverse_shift(P, in, dt); + + return out.lpz; +} + +static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + PJ_COORD out = in; + double dt; + + + if (Q->t_obs != HUGE_VAL) { + dt = Q->t_epoch - Q->t_obs; + } else { + dt = Q->t_epoch - in.xyzt.t; + } + + out.xyz = reverse_shift(P, in.xyz, dt); + return out; +} + +static PJ *destructor(PJ *P, int errlev) { + if (nullptr==P) + return nullptr; + + if (nullptr==P->opaque) + return pj_default_destructor (P, errlev); + + if (static_cast(P->opaque)->cart) + static_cast(P->opaque)->cart->destructor (static_cast(P->opaque)->cart, errlev); + + return pj_default_destructor(P, errlev); +} + + +PJ *TRANSFORMATION(deformation,1) { + int has_xy_grids = 0; + int has_z_grids = 0; + struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return destructor(P, ENOMEM); + P->opaque = (void *) Q; + + Q->cart = proj_create(P->ctx, "+proj=cart"); + if (Q->cart == nullptr) + return destructor(P, ENOMEM); + + /* inherit ellipsoid definition from P to Q->cart */ + pj_inherit_ellipsoid_def (P, Q->cart); + + has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i; + has_z_grids = pj_param(P->ctx, P->params, "tz_grids").i; + + /* Build gridlists. Both horizontal and vertical grids are mandatory. */ + if (!has_xy_grids || !has_z_grids) { + proj_log_error(P, "deformation: Both +xy_grids and +z_grids should be specified."); + return destructor(P, PJD_ERR_NO_ARGS ); + } + + proj_hgrid_init(P, "xy_grids"); + if (proj_errno(P)) { + proj_log_error(P, "deformation: could not find requested xy_grid(s)."); + return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + proj_vgrid_init(P, "z_grids"); + if (proj_errno(P)) { + proj_log_error(P, "deformation: could not find requested z_grid(s)."); + return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + Q->t_obs = HUGE_VAL; + if (pj_param(P->ctx, P->params, "tt_obs").i) { + Q->t_obs = pj_param(P->ctx, P->params, "dt_obs").f; + } + + if (pj_param(P->ctx, P->params, "tt_epoch").i) { + Q->t_epoch = pj_param(P->ctx, P->params, "dt_epoch").f; + } else { + proj_log_error(P, "deformation: +t_epoch parameter missing."); + return destructor(P, PJD_ERR_MISSING_ARGS); + } + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = nullptr; + P->inv = nullptr; + + P->left = PJ_IO_UNITS_CARTESIAN; + P->right = PJ_IO_UNITS_CARTESIAN; + P->destructor = destructor; + + return P; +} + diff --git a/src/transformations/helmert.cpp b/src/transformations/helmert.cpp new file mode 100644 index 00000000..4a3abf4e --- /dev/null +++ b/src/transformations/helmert.cpp @@ -0,0 +1,755 @@ +/*********************************************************************** + + 3-, 4-and 7-parameter shifts, and their 6-, 8- + and 14-parameter kinematic counterparts. + + Thomas Knudsen, 2016-05-24 + +************************************************************************ + + 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. + +************************************************************************ + +Thomas Knudsen, thokn@sdfe.dk, 2016-05-24/06-05 +Kristian Evers, kreve@sdfe.dk, 2017-05-01 +Even Rouault, even.roault@spatialys.com +Last update: 2018-10-26 + +************************************************************************ +* Copyright (c) 2016, Thomas Knudsen / SDFE +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +* DEALINGS IN THE SOFTWARE. +* +***********************************************************************/ + +#define PJ_LIB__ + +#include +#include + +#include "proj_internal.h" +#include "projects.h" +#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); + + + +/***********************************************************************/ +namespace { // anonymous namespace +struct pj_opaque_helmert { +/************************************************************************ + Projection specific elements for the "helmert" PJ object +************************************************************************/ + XYZ xyz; + XYZ xyz_0; + XYZ dxyz; + XYZ refp; + PJ_OPK opk; + PJ_OPK opk_0; + PJ_OPK dopk; + double scale; + double scale_0; + double dscale; + double theta; + double theta_0; + double dtheta; + double R[3][3]; + double t_epoch, t_obs; + int no_rotation, exact, fourparam; + int is_position_vector; /* 1 = position_vector, 0 = coordinate_frame */ +}; +} // anonymous namespace + + +/* Make the maths of the rotation operations somewhat more readable and textbook like */ +#define R00 (Q->R[0][0]) +#define R01 (Q->R[0][1]) +#define R02 (Q->R[0][2]) + +#define R10 (Q->R[1][0]) +#define R11 (Q->R[1][1]) +#define R12 (Q->R[1][2]) + +#define R20 (Q->R[2][0]) +#define R21 (Q->R[2][1]) +#define R22 (Q->R[2][2]) + +/**************************************************************************/ +static void update_parameters(PJ *P) { +/*************************************************************************** + + Update transformation parameters. + --------------------------------- + + The 14-parameter Helmert transformation is at it's core the same as the + 7-parameter transformation, since the transformation parameters are + projected forward or backwards in time via the rate of changes of the + parameters. The transformation parameters are calculated for a specific + epoch before the actual Helmert transformation is carried out. + + The transformation parameters are updated with the following equation [0]: + + . + P(t) = P(EPOCH) + P * (t - EPOCH) + + . + where EPOCH is the epoch indicated in the above table and P is the rate + of that parameter. + + [0] http://itrf.ign.fr/doc_ITRF/Transfo-ITRF2008_ITRFs.txt + +*******************************************************************************/ + + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + double dt = Q->t_obs - Q->t_epoch; + + Q->xyz.x = Q->xyz_0.x + Q->dxyz.x * dt; + Q->xyz.y = Q->xyz_0.y + Q->dxyz.y * dt; + Q->xyz.z = Q->xyz_0.z + Q->dxyz.z * dt; + + Q->opk.o = Q->opk_0.o + Q->dopk.o * dt; + Q->opk.p = Q->opk_0.p + Q->dopk.p * dt; + Q->opk.k = Q->opk_0.k + Q->dopk.k * dt; + + Q->scale = Q->scale_0 + Q->dscale * dt; + + Q->theta = Q->theta_0 + Q->dtheta * dt; + + /* debugging output */ + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) { + proj_log_trace(P, "Transformation parameters for observation " + "t_obs=%g (t_epoch=%g):", Q->t_obs, Q->t_epoch); + proj_log_trace(P, "x: %g", Q->xyz.x); + proj_log_trace(P, "y: %g", Q->xyz.y); + proj_log_trace(P, "z: %g", Q->xyz.z); + proj_log_trace(P, "s: %g", Q->scale*1e-6); + proj_log_trace(P, "rx: %g", Q->opk.o); + proj_log_trace(P, "ry: %g", Q->opk.p); + proj_log_trace(P, "rz: %g", Q->opk.k); + proj_log_trace(P, "theta: %g", Q->theta); + } +} + +/**************************************************************************/ +static void build_rot_matrix(PJ *P) { +/*************************************************************************** + + Build rotation matrix. + ---------------------- + + Here we rename rotation indices from omega, phi, kappa (opk), to + fi (i.e. phi), theta, psi (ftp), in order to reduce the mental agility + needed to implement the expression for the rotation matrix derived over + at https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions + The relevant section is Euler angles ( z-’-x" intrinsic) -> Rotation matrix + + By default small angle approximations are used: + The matrix elements are approximated by expanding the trigonometric + functions to linear order (i.e. cos(x) = 1, sin(x) = x), and discarding + products of second order. + + This was a useful hack when calculating by hand was the only option, + but in general, today, should be avoided because: + + 1. It does not save much computation time, as the rotation matrix + is built only once and probably used many times (except when + transforming spatio-temporal coordinates). + + 2. The error induced may be too large for ultra high accuracy + applications: the Earth is huge and the linear error is + approximately the angular error multiplied by the Earth radius. + + However, in many cases the approximation is necessary, since it has + been used historically: Rotation angles from older published datum + shifts may actually be a least squares fit to the linearized rotation + approximation, hence not being strictly valid for deriving the exact + rotation matrix. In fact, most publicly available transformation + parameters are based on the approximate Helmert transform, which is why + we use that as the default setting, even though it is more correct to + use the exact form of the equations. + + So in order to fit historically derived coordinates, the access to + the approximate rotation matrix is necessary - at least in principle. + + Also, when using any published datum transformation information, one + should always check which convention (exact or approximate rotation + matrix) is expected, and whether the induced error for selecting + the opposite convention is acceptable (which it often is). + + + Sign conventions + ---------------- + + Take care: Two different sign conventions exist for the rotation terms. + + Conceptually they relate to whether we rotate the coordinate system + or the "position vector" (the vector going from the coordinate system + origin to the point being transformed, i.e. the point coordinates + interpreted as vector coordinates). + + Switching between the "position vector" and "coordinate system" + conventions is simply a matter of switching the sign of the rotation + angles, which algebraically also translates into a transposition of + the rotation matrix. + + Hence, as geodetic constants should preferably be referred to exactly + as published, the "convention" option provides the ability to switch + between the conventions. + +***************************************************************************/ + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + + double f, t, p; /* phi/fi , theta, psi */ + double cf, ct, cp; /* cos (fi, theta, psi) */ + double sf, st, sp; /* sin (fi, theta, psi) */ + + /* rename (omega, phi, kappa) to (fi, theta, psi) */ + f = Q->opk.o; + t = Q->opk.p; + p = Q->opk.k; + + /* Those equations are given assuming coordinate frame convention. */ + /* For the position vector convention, we transpose the matrix just after. */ + if (Q->exact) { + cf = cos(f); + sf = sin(f); + ct = cos(t); + st = sin(t); + cp = cos(p); + sp = sin(p); + + + R00 = ct*cp; + R01 = cf*sp + sf*st*cp; + R02 = sf*sp - cf*st*cp; + + R10 = -ct*sp; + R11 = cf*cp - sf*st*sp; + R12 = sf*cp + cf*st*sp; + + R20 = st; + R21 = -sf*ct; + R22 = cf*ct; + } else{ + R00 = 1; + R01 = p; + R02 = -t; + + R10 = -p; + R11 = 1; + R12 = f; + + R20 = t; + R21 = -f; + R22 = 1; + } + + + /* + For comparison: Description from Engsager/Poder implementation + in set_dtm_1.c (trlib) + + DATUM SHIFT: + TO = scale * ROTZ * ROTY * ROTX * FROM + TRANSLA + + ( cz sz 0) (cy 0 -sy) (1 0 0) + ROTZ=(-sz cz 0), ROTY=(0 1 0), ROTX=(0 cx sx) + ( 0 0 1) (sy 0 cy) (0 -sx cx) + + trp->r11 = cos_ry*cos_rz; + trp->r12 = cos_rx*sin_rz + sin_rx*sin_ry*cos_rz; + trp->r13 = sin_rx*sin_rz - cos_rx*sin_ry*cos_rz; + + trp->r21 = -cos_ry*sin_rz; + trp->r22 = cos_rx*cos_rz - sin_rx*sin_ry*sin_rz; + trp->r23 = sin_rx*cos_rz + cos_rx*sin_ry*sin_rz; + + trp->r31 = sin_ry; + trp->r32 = -sin_rx*cos_ry; + trp->r33 = cos_rx*cos_ry; + + trp->scale = 1.0 + scale; + */ + + + if (Q->is_position_vector) { + double r; + r = R01; R01 = R10; R10 = r; + r = R02; R02 = R20; R20 = r; + r = R12; R12 = R21; R21 = r; + } + + /* some debugging output */ + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) { + proj_log_trace(P, "Rotation Matrix:"); + proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R00, R01, R02); + proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R10, R11, R12); + proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R20, R21, R22); + } +} + + + + +/***********************************************************************/ +static XY helmert_forward (LP lp, PJ *P) { +/***********************************************************************/ + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + double x, y, cr, sr; + point.lp = lp; + + cr = cos(Q->theta) * Q->scale; + sr = sin(Q->theta) * Q->scale; + x = point.xy.x; + y = point.xy.y; + + point.xy.x = cr*x + sr*y + Q->xyz_0.x; + point.xy.y = -sr*x + cr*y + Q->xyz_0.y; + + return point.xy; +} + + +/***********************************************************************/ +static LP helmert_reverse (XY xy, PJ *P) { +/***********************************************************************/ + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + double x, y, sr, cr; + point.xy = xy; + + cr = cos(Q->theta) / Q->scale; + sr = sin(Q->theta) / Q->scale; + x = point.xy.x - Q->xyz_0.x; + y = point.xy.y - Q->xyz_0.y; + + point.xy.x = x*cr - y*sr; + point.xy.y = x*sr + y*cr; + + return point.lp; +} + + +/***********************************************************************/ +static XYZ helmert_forward_3d (LPZ lpz, PJ *P) { +/***********************************************************************/ + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + double X, Y, Z, scale; + + point.lpz = lpz; + + if (Q->fourparam) { + point.xy = helmert_forward(point.lp, P); + return point.xyz; + } + + if (Q->no_rotation) { + point.xyz.x = lpz.lam + Q->xyz.x; + point.xyz.y = lpz.phi + Q->xyz.y; + point.xyz.z = lpz.z + Q->xyz.z; + return point.xyz; + } + + scale = 1 + Q->scale * 1e-6; + + 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; /* for Molodensky-Badekas, Q->xyz already incorporates the Q->refp offset */ + point.xyz.y += Q->xyz.y; + point.xyz.z += Q->xyz.z; + + return point.xyz; +} + + +/***********************************************************************/ +static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) { +/***********************************************************************/ + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + double X, Y, Z, scale; + + point.xyz = xyz; + + if (Q->fourparam) { + point.lp = helmert_reverse(point.xy, P); + return point.lpz; + } + + if (Q->no_rotation) { + point.xyz.x = xyz.x - Q->xyz.x; + point.xyz.y = xyz.y - Q->xyz.y; + point.xyz.z = xyz.z - Q->xyz.z; + return point.lpz; + } + + scale = 1 + Q->scale * 1e-6; + + /* Unscale and deoffset */ + X = (xyz.x - Q->xyz.x) / scale; + Y = (xyz.y - Q->xyz.y) / scale; + Z = (xyz.z - Q->xyz.z) / scale; + + /* Inverse rotation through transpose multiplication */ + 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; +} + + +static PJ_COORD helmert_forward_4d (PJ_COORD point, PJ *P) { + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + + /* We only need to rebuild the rotation matrix if the + * observation time is different from the last call */ + double t_obs = (point.xyzt.t == HUGE_VAL) ? Q->t_epoch : point.xyzt.t; + if (t_obs != Q->t_obs) { + Q->t_obs = t_obs; + update_parameters(P); + build_rot_matrix(P); + } + + point.xyz = helmert_forward_3d (point.lpz, P); + + return point; +} + + +static PJ_COORD helmert_reverse_4d (PJ_COORD point, PJ *P) { + struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; + + /* We only need to rebuild the rotation matrix if the + * observation time is different from the last call */ + double t_obs = (point.xyzt.t == HUGE_VAL) ? Q->t_epoch : point.xyzt.t; + if (t_obs != Q->t_obs) { + Q->t_obs = t_obs; + update_parameters(P); + build_rot_matrix(P); + } + + point.lpz = helmert_reverse_3d (point.xyz, P); + + return point; +} + +/* Arcsecond to radians */ +#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) + + +static PJ* init_helmert_six_parameters(PJ* P) { + struct pj_opaque_helmert *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_helmert))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = (void *) Q; + + /* In most cases, we work on 3D cartesian coordinates */ + P->left = PJ_IO_UNITS_CARTESIAN; + P->right = PJ_IO_UNITS_CARTESIAN; + + /* 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 nullptr; + } + + /* 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. " + "Use convention=position_vector/coordinate_frame"); + return pj_default_destructor (P, PJD_ERR_INVALID_ARG); + } + + /* Support the classic PROJ towgs84 parameter, but allow later overrides.*/ + /* Note that if towgs84 is specified, the datum_params array is set up */ + /* for us automagically by the pj_datum_set call in pj_init_ctx */ + if (pj_param_exists (P->params, "towgs84")) { + Q->xyz_0.x = P->datum_params[0]; + Q->xyz_0.y = P->datum_params[1]; + Q->xyz_0.z = P->datum_params[2]; + + Q->opk_0.o = P->datum_params[3]; + Q->opk_0.p = P->datum_params[4]; + Q->opk_0.k = P->datum_params[5]; + + /* We must undo conversion to absolute scale from pj_datum_set */ + if (0==P->datum_params[6]) + Q->scale_0 = 0; + else + Q->scale_0 = (P->datum_params[6] - 1) * 1e6; + } + + 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; + Q->scale_0 = 1.0; /* default scale for the 4-param shift */ + } + + /* Scale */ + if (pj_param (P->ctx, P->params, "ts").i) { + Q->scale_0 = pj_param (P->ctx, P->params, "ds").f; + if (pj_param (P->ctx, P->params, "ttheta").i && Q->scale_0 == 0.0) + return pj_default_destructor (P, PJD_ERR_INVALID_SCALE); + } + + /* Translation rates */ + if (pj_param(P->ctx, P->params, "tdx").i) + Q->dxyz.x = pj_param (P->ctx, P->params, "ddx").f; + + if (pj_param(P->ctx, P->params, "tdy").i) + Q->dxyz.y = pj_param (P->ctx, P->params, "ddy").f; + + if (pj_param(P->ctx, P->params, "tdz").i) + Q->dxyz.z = pj_param (P->ctx, P->params, "ddz").f; + + /* Rotations rates */ + if (pj_param (P->ctx, P->params, "tdrx").i) + Q->dopk.o = pj_param (P->ctx, P->params, "ddrx").f * ARCSEC_TO_RAD; + + if (pj_param (P->ctx, P->params, "tdry").i) + Q->dopk.p = pj_param (P->ctx, P->params, "ddry").f * ARCSEC_TO_RAD; + + if (pj_param (P->ctx, P->params, "tdrz").i) + Q->dopk.k = pj_param (P->ctx, P->params, "ddrz").f * ARCSEC_TO_RAD; + + if (pj_param (P->ctx, P->params, "tdtheta").i) + Q->dtheta = pj_param (P->ctx, P->params, "ddtheta").f * ARCSEC_TO_RAD; + + /* Scale rate */ + if (pj_param (P->ctx, P->params, "tds").i) + Q->dscale = pj_param (P->ctx, P->params, "dds").f; + + + /* Epoch */ + if (pj_param(P->ctx, P->params, "tt_epoch").i) + Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f; + + if (pj_param(P->ctx, P->params, "tt_obs").i) + Q->t_obs = pj_param (P->ctx, P->params, "dt_obs").f; + + Q->xyz = Q->xyz_0; + Q->opk = Q->opk_0; + Q->scale = Q->scale_0; + Q->theta = Q->theta_0; + + if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0) && + (Q->dopk.o==0) && (Q->dopk.p==0) && (Q->dopk.k==0)) { + Q->no_rotation = 1; + } + + if( !read_convention(P) ) { + return nullptr; + } + + /* Let's help with debugging */ + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) { + proj_log_debug(P, "Helmert parameters:"); + proj_log_debug(P, "x= %8.5f y= %8.5f z= %8.5f", Q->xyz.x, Q->xyz.y, Q->xyz.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->no_rotation ? "" : + Q->is_position_vector ? " convention=position_vector" : + " convention=coordinate_frame"); + proj_log_debug(P, "dx= %8.5f dy= %8.5f dz= %8.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z); + proj_log_debug(P, "drx=%8.5f dry=%8.5f drz=%8.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k); + proj_log_debug(P, "ds= %8.5f t_epoch=%8.5f t_obs=%8.5f", Q->dscale, Q->t_epoch, Q->t_obs); + } + + if (Q->no_rotation) { + return P; + } + + update_parameters(P); + build_rot_matrix(P); + + return P; +} + + +/***********************************************************************/ +PJ *TRANSFORMATION(molobadekas, 0) { +/***********************************************************************/ + + struct pj_opaque_helmert *Q; + + if( !init_helmert_six_parameters(P) ) { + return nullptr; + } + + 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 nullptr; + } + + /* 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/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp new file mode 100644 index 00000000..f0e57251 --- /dev/null +++ b/src/transformations/hgridshift.cpp @@ -0,0 +1,133 @@ +#define PJ_LIB__ + +#include +#include +#include +#include + +#include "proj_internal.h" +#include "projects.h" + +PROJ_HEAD(hgridshift, "Horizontal grid shift"); + +namespace { // anonymous namespace +struct pj_opaque_hgridshift { + double t_final; + double t_epoch; +}; +} // anonymous namespace + +static XYZ forward_3d(LPZ lpz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lpz = lpz; + + if (P->gridlist != nullptr) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD); + } + + return point.xyz; +} + + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xyz = xyz; + + if (P->gridlist != nullptr) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + point.lp = proj_hgrid_apply(P, point.lp, PJ_INV); + } + + return point.lpz; +} + +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { + struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque; + PJ_COORD point = obs; + + /* If transformation is not time restricted, we always call it */ + if (Q->t_final==0 || Q->t_epoch==0) { + point.xyz = forward_3d (obs.lpz, P); + return point; + } + + /* Time restricted - only apply transform if within time bracket */ + if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) + point.xyz = forward_3d (obs.lpz, P); + + + return point; +} + +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { + struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque; + PJ_COORD point = obs; + + /* If transformation is not time restricted, we always call it */ + if (Q->t_final==0 || Q->t_epoch==0) { + point.lpz = reverse_3d (obs.xyz, P); + return point; + } + + /* Time restricted - only apply transform if within time bracket */ + if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) + point.lpz = reverse_3d (obs.xyz, P); + + return point; +} + + +PJ *TRANSFORMATION(hgridshift,0) { + struct pj_opaque_hgridshift *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_hgridshift))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = nullptr; + P->inv = nullptr; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + if (0==pj_param(P->ctx, P->params, "tgrids").i) { + proj_log_error(P, "hgridshift: +grids parameter missing."); + return pj_default_destructor (P, PJD_ERR_NO_ARGS); + } + + /* TODO: Refactor into shared function that can be used */ + /* by both vgridshift and hgridshift */ + if (pj_param(P->ctx, P->params, "tt_final").i) { + Q->t_final = pj_param (P->ctx, P->params, "dt_final").f; + if (Q->t_final == 0) { + /* a number wasn't passed to +t_final, let's see if it was "now" */ + /* and set the time accordingly. */ + if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) { + time_t now; + struct tm *date; + time(&now); + date = localtime(&now); + Q->t_final = 1900.0 + date->tm_year + date->tm_yday/365.0; + } + } + } + + if (pj_param(P->ctx, P->params, "tt_epoch").i) + Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f; + + + proj_hgrid_init(P, "grids"); + /* Was gridlist compiled properly? */ + if ( proj_errno(P) ) { + proj_log_error(P, "hgridshift: could not find required grid(s)."); + return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + return P; +} diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp new file mode 100644 index 00000000..73977de6 --- /dev/null +++ b/src/transformations/horner.cpp @@ -0,0 +1,513 @@ +/*********************************************************************** + + Interfacing to a classic piece of geodetic software + +************************************************************************ + + gen_pol is a highly efficient, classic implementation of a generic + 2D Horner's Scheme polynomial evaluation routine by Knud Poder and + Karsten Engsager, originating in the vivid geodetic environment at + what was then (1960-ish) the Danish Geodetic Institute. + + The original Poder/Engsager gen_pol implementation (where + the polynomial degree and two sets of polynomial coefficients + are packed together in one compound array, handled via a plain + double pointer) is compelling and "true to the code history": + + It has a beautiful classical 1960s ring to it, not unlike the + original fft implementations, which revolutionized spectral + analysis in twenty lines of code. + + The Poder coding sound, as classic 1960s as Phil Spector's Wall + of Sound, is beautiful and inimitable. + + On the other hand: For the uninitiated, the gen_pol code is hard + to follow, despite being compact. + + Also, since adding metadata and improving maintainability + of the code are among the implied goals of a current SDFE/DTU Space + project, the material in this file introduces a version with a + more modern (or at least 1990s) look, introducing a "double 2D + polynomial" data type, HORNER. + + Despite introducing a new data type for handling the polynomial + coefficients, great care has been taken to keep the coefficient + array organization identical to that of gen_pol. + + Hence, on one hand, the HORNER data type helps improving the + long term maintainability of the code by making the data + organization more mentally accessible. + + On the other hand, it allows us to preserve the business end of + the original gen_pol implementation - although not including the + famous "Poder dual autocheck" in all its enigmatic elegance. + + ********************************************************************** + + The material included here was written by Knud Poder, starting + around 1960, and Karsten Engsager, starting around 1970. It was + originally written in Algol 60, later (1980s) reimplemented in C. + + The HORNER data type interface, and the organization as a header + library was implemented by Thomas Knudsen, starting around 2015. + + *********************************************************************** + * + * Copyright (c) 2016, SDFE http://www.sdfe.dk / Thomas Knudsen / Karsten Engsager + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + * + *****************************************************************************/ + +#define PJ_LIB__ + +#include +#include +#include +#include +#include + +#include "proj_internal.h" +#include "projects.h" + +PROJ_HEAD(horner, "Horner polynomial evaluation"); + +/* make horner.h interface with proj's memory management */ +#define horner_dealloc(x) pj_dealloc(x) +#define horner_calloc(n,x) pj_calloc(n,x) + +namespace { // anonymous namespace +struct horner { + int uneg; /* u axis negated? */ + int vneg; /* v axis negated? */ + int order; /* maximum degree of polynomium */ + int coefs; /* number of coefficients for each polynomium */ + double range; /* radius of the region of validity */ + + double *fwd_u; /* coefficients for the forward transformations */ + double *fwd_v; /* i.e. latitude/longitude to northing/easting */ + + double *inv_u; /* coefficients for the inverse transformations */ + double *inv_v; /* i.e. northing/easting to latitude/longitude */ + + double *fwd_c; /* coefficients for the complex forward transformations */ + double *inv_c; /* coefficients for the complex inverse transformations */ + + UV *fwd_origin; /* False longitude/latitude */ + UV *inv_origin; /* False easting/northing */ +}; +} // anonymous namespace + +typedef struct horner HORNER; +static UV horner_func (const HORNER *transformation, PJ_DIRECTION direction, UV position); +static HORNER *horner_alloc (size_t order, int complex_polynomia); +static void horner_free (HORNER *h); + +/* e.g. degree = 2: a + bx + cy + dxx + eyy + fxy, i.e. 6 coefficients */ +#define horner_number_of_coefficients(order) \ + (((order + 1)*(order + 2)/2)) + + +static void horner_free (HORNER *h) { + horner_dealloc (h->inv_v); + horner_dealloc (h->inv_u); + horner_dealloc (h->fwd_v); + horner_dealloc (h->fwd_u); + horner_dealloc (h->fwd_c); + horner_dealloc (h->inv_c); + horner_dealloc (h->fwd_origin); + horner_dealloc (h->inv_origin); + horner_dealloc (h); +} + + +static HORNER *horner_alloc (size_t order, int complex_polynomia) { + /* size_t is unsigned, so we need not check for order > 0 */ + int n = (int)horner_number_of_coefficients(order); + int polynomia_ok = 0; + HORNER *h = static_cast(horner_calloc (1, sizeof (HORNER))); + + if (nullptr==h) + return nullptr; + + if (complex_polynomia) + n = 2*(int)order + 2; + h->order = (int)order; + h->coefs = n; + + if (complex_polynomia) { + h->fwd_c = static_cast(horner_calloc (n, sizeof(double))); + h->inv_c = static_cast(horner_calloc (n, sizeof(double))); + if (h->fwd_c && h->inv_c) + polynomia_ok = 1; + } + else { + h->fwd_u = static_cast(horner_calloc (n, sizeof(double))); + h->fwd_v = static_cast(horner_calloc (n, sizeof(double))); + h->inv_u = static_cast(horner_calloc (n, sizeof(double))); + h->inv_v = static_cast(horner_calloc (n, sizeof(double))); + if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v) + polynomia_ok = 1; + } + + h->fwd_origin = static_cast(horner_calloc (1, sizeof(UV))); + h->inv_origin = static_cast(horner_calloc (1, sizeof(UV))); + + if (polynomia_ok && h->fwd_origin && h->inv_origin) + return h; + + /* safe, since all pointers are null-initialized (by calloc) */ + horner_free (h); + return nullptr; +} + + + + +/**********************************************************************/ +static UV horner_func (const HORNER *transformation, PJ_DIRECTION direction, UV position) { +/*********************************************************************** + +A reimplementation of the classic Engsager/Poder 2D Horner polynomial +evaluation engine "gen_pol". + +This version omits the inimitable Poder "dual autocheck"-machinery, +which here is intended to be implemented at a higher level of the +library: We separate the polynomial evaluation from the quality +control (which, given the limited MTBF for "computing machinery", +typical when Knud Poder invented the dual autocheck method, +was not defensible at that time). + +Another difference from the original version is that we return the +result on the stack, rather than accepting pointers to result variables +as input. This results in code that is easy to read: + + projected = horner (s34j, 1, geographic); + geographic = horner (s34j, -1, projected ); + +and experiments have shown that on contemporary architectures, the time +taken for returning even comparatively large objects on the stack (and +the UV is not that large - typically only 16 bytes) is negligibly +different from passing two pointers (i.e. typically also 16 bytes) the +other way. + +The polynomium has the form: + +P = sum (i = [0 : order]) + sum (j = [0 : order - i]) + pow(par_1, i) * pow(par_2, j) * coef(index(order, i, j)) + +For numerical stability, the summation is carried out backwards, +summing the tiny high order elements first. + +***********************************************************************/ + + /* These variable names follow the Engsager/Poder implementation */ + int sz; /* Number of coefficients per polynomial */ + double *tcx, *tcy; /* Coefficient pointers */ + double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ + double n, e; + UV uv_error; + uv_error.u = uv_error.v = HUGE_VAL; + + if (nullptr==transformation) + return uv_error; + + /* Check for valid value of direction (-1, 0, 1) */ + switch (direction) { + case PJ_IDENT: /* no-op */ + return position; + case PJ_FWD: /* forward */ + case PJ_INV: /* inverse */ + break; + default: /* invalid */ + errno = EINVAL; + return uv_error; + } + + /* Prepare for double Horner */ + sz = horner_number_of_coefficients(transformation->order); + range = transformation->range; + + + if (direction==PJ_FWD) { /* forward */ + tcx = transformation->fwd_u + sz; + tcy = transformation->fwd_v + sz; + e = position.u - transformation->fwd_origin->u; + n = position.v - transformation->fwd_origin->v; + } else { /* inverse */ + tcx = transformation->inv_u + sz; + tcy = transformation->inv_v + sz; + e = position.u - transformation->inv_origin->u; + n = position.v - transformation->inv_origin->v; + } + + if ((fabs(n) > range) || (fabs(e) > range)) { + errno = EDOM; + return uv_error; + } + + /* The melody of this block is straight out of the great Engsager/Poder songbook */ + else { + int g = transformation->order; + int r = g, c; + double u, v, N, E; + + /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */ + N = *--tcy; + E = *--tcx; + for (; r > 0; r--) { + u = *--tcy; + v = *--tcx; + for (c = g; c >= r; c--) { + u = n*u + *--tcy; + v = e*v + *--tcx; + } + N = e*N + u; + E = n*E + v; + } + + position.u = E; + position.v = N; + } + + return position; +} + + + + + + + +static PJ_COORD horner_forward_4d (PJ_COORD point, PJ *P) { + point.uv = horner_func ((HORNER *) P->opaque, PJ_FWD, point.uv); + return point; +} + +static PJ_COORD horner_reverse_4d (PJ_COORD point, PJ *P) { + point.uv = horner_func ((HORNER *) P->opaque, PJ_INV, point.uv); + return point; +} + + + + +/**********************************************************************/ +static UV complex_horner (const HORNER *transformation, PJ_DIRECTION direction, UV position) { +/*********************************************************************** + +A reimplementation of a classic Engsager/Poder Horner complex +polynomial evaluation engine. + +***********************************************************************/ + + /* These variable names follow the Engsager/Poder implementation */ + int sz; /* Number of coefficients */ + double *c, *cb; /* Coefficient pointers */ + double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ + double n, e, w, N, E; + UV uv_error; + uv_error.u = uv_error.v = HUGE_VAL; + + if (nullptr==transformation) + return uv_error; + + /* Check for valid value of direction (-1, 0, 1) */ + switch (direction) { + case PJ_IDENT: /* no-op */ + return position; + case PJ_FWD: /* forward */ + case PJ_INV: /* inverse */ + break; + default: /* invalid */ + errno = EINVAL; + return uv_error; + } + + /* Prepare for double Horner */ + sz = 2*transformation->order + 2; + range = transformation->range; + + if (direction==PJ_FWD) { /* forward */ + cb = transformation->fwd_c; + c = cb + sz; + e = position.u - transformation->fwd_origin->u; + n = position.v - transformation->fwd_origin->v; + if (transformation->uneg) + e = -e; + if (transformation->vneg) + n = -n; + } else { /* inverse */ + cb = transformation->inv_c; + c = cb + sz; + e = position.u - transformation->inv_origin->u; + n = position.v - transformation->inv_origin->v; + if (transformation->uneg) + e = -e; + if (transformation->vneg) + n = -n; + } + + if ((fabs(n) > range) || (fabs(e) > range)) { + errno = EDOM; + return uv_error; + } + + /* Everything's set up properly - now do the actual polynomium evaluation */ + E = *--c; + N = *--c; + while (c > cb) { + w = n*E + e*N + *--c; + N = n*N - e*E + *--c; + E = w; + } + + position.u = E; + position.v = N; + return position; +} + + + +static PJ_COORD complex_horner_forward_4d (PJ_COORD point, PJ *P) { + point.uv = complex_horner ((HORNER *) P->opaque, PJ_FWD, point.uv); + return point; +} + +static PJ_COORD complex_horner_reverse_4d (PJ_COORD point, PJ *P) { + point.uv = complex_horner ((HORNER *) P->opaque, PJ_INV, point.uv); + return point; +} + + +static PJ *horner_freeup (PJ *P, int errlev) { /* Destructor */ + if (nullptr==P) + return nullptr; + if (nullptr==P->opaque) + return pj_default_destructor (P, errlev); + horner_free ((HORNER *) P->opaque); + P->opaque = nullptr; + return pj_default_destructor (P, errlev); +} + + +static int parse_coefs (PJ *P, double *coefs, const char *param, int ncoefs) { + char *buf, *init, *next = nullptr; + int i; + + buf = static_cast(pj_calloc (strlen (param) + 2, sizeof(char))); + if (nullptr==buf) { + proj_log_error (P, "Horner: No memory left"); + return 0; + } + + sprintf (buf, "t%s", param); + if (0==pj_param (P->ctx, P->params, buf).i) { + pj_dealloc (buf); + return 0; + } + sprintf (buf, "s%s", param); + init = pj_param(P->ctx, P->params, buf).s; + pj_dealloc (buf); + + for (i = 0; i < ncoefs; i++) { + if (i > 0) { + if ( next == nullptr || ','!=*next) { + proj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs); + return 0; + } + init = ++next; + } + coefs[i] = pj_strtod (init, &next); + } + return 1; +} + + +/*********************************************************************/ +PJ *PROJECTION(horner) { +/*********************************************************************/ + int degree = 0, n, complex_polynomia = 0; + HORNER *Q; + P->fwd4d = horner_forward_4d; + P->inv4d = horner_reverse_4d; + P->fwd3d = nullptr; + P->inv3d = nullptr; + P->fwd = nullptr; + P->inv = nullptr; + P->left = P->right = PJ_IO_UNITS_PROJECTED; + P->destructor = horner_freeup; + + /* Polynomial degree specified? */ + if (pj_param (P->ctx, P->params, "tdeg").i) { /* degree specified? */ + degree = pj_param(P->ctx, P->params, "ideg").i; + if (degree < 0 || degree > 10000) { + /* What are reasonable minimum and maximums for degree? */ + proj_log_debug (P, "Horner: Degree is unreasonable: %d", degree); + return horner_freeup (P, PJD_ERR_INVALID_ARG); + } + } else { + proj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)"); + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + } + + if (pj_param (P->ctx, P->params, "tfwd_c").i || pj_param (P->ctx, P->params, "tinv_c").i) /* complex polynomium? */ + complex_polynomia = 1; + + Q = horner_alloc (degree, complex_polynomia); + if (Q == nullptr) + return horner_freeup (P, ENOMEM); + P->opaque = Q; + + if (complex_polynomia) { + /* Westings and/or southings? */ + Q->uneg = pj_param_exists (P->params, "uneg") ? 1 : 0; + Q->vneg = pj_param_exists (P->params, "vneg") ? 1 : 0; + + n = 2*degree + 2; + if (0==parse_coefs (P, Q->fwd_c, "fwd_c", n)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + if (0==parse_coefs (P, Q->inv_c, "inv_c", n)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + P->fwd4d = complex_horner_forward_4d; + P->inv4d = complex_horner_reverse_4d; + } + + else { + n = horner_number_of_coefficients (degree); + if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + if (0==parse_coefs (P, Q->inv_u, "inv_u", n)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + if (0==parse_coefs (P, Q->inv_v, "inv_v", n)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + } + + if (0==parse_coefs (P, (double *)(Q->fwd_origin), "fwd_origin", 2)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) + return horner_freeup (P, PJD_ERR_MISSING_ARGS); + if (0==parse_coefs (P, &Q->range, "range", 1)) + Q->range = 500000; + + return P; +} diff --git a/src/transformations/molodensky.cpp b/src/transformations/molodensky.cpp new file mode 100644 index 00000000..91743fda --- /dev/null +++ b/src/transformations/molodensky.cpp @@ -0,0 +1,329 @@ +/*********************************************************************** + + (Abridged) Molodensky Transform + + Kristian Evers, 2017-07-07 + +************************************************************************ + + Implements the (abridged) Molodensky transformations for 2D and 3D + data. + + Primarily useful for implementation of datum shifts in transformation + pipelines. + + The code in this file is mostly based on + + The Standard and Abridged Molodensky Coordinate Transformation + Formulae, 2004, R.E. Deaking, + http://www.mygeodesy.id.au/documents/Molodensky%20V2.pdf + + + +************************************************************************ +* Copyright (c) 2017, Kristian Evers / SDFE +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +* DEALINGS IN THE SOFTWARE. +* +***********************************************************************/ +#define PJ_LIB__ + +#include +#include + +#include "proj.h" +#include "proj_internal.h" +#include "projects.h" + +PROJ_HEAD(molodensky, "Molodensky transform"); + +static XYZ forward_3d(LPZ lpz, PJ *P); +static LPZ reverse_3d(XYZ xyz, PJ *P); + +namespace { // anonymous namespace +struct pj_opaque_molodensky { + double dx; + double dy; + double dz; + double da; + double df; + int abridged; +}; +} // anonymous namespace + + +static double RN (double a, double es, double phi) { +/********************************************************** + N(phi) - prime vertical radius of curvature + ------------------------------------------- + + This is basically the same function as in PJ_cart.c + should probably be refactored into it's own file at some + point. + +**********************************************************/ + double s = sin(phi); + if (es==0) + return a; + + return a / sqrt (1 - es*s*s); +} + + +static double RM (double a, double es, double phi) { +/********************************************************** + M(phi) - Meridian radius of curvature + ------------------------------------- + + Source: + + E.J Krakiwsky & D.B. Thomson, 1974, + GEODETIC POSITION COMPUTATIONS, + + Fredericton NB, Canada: + University of New Brunswick, + Department of Geodesy and Geomatics Engineering, + Lecture Notes No. 39, + 99 pp. + + http://www2.unb.ca/gge/Pubs/LN39.pdf + +**********************************************************/ + double s = sin(phi); + if (es==0) + return a; + + /* eq. 13a */ + if (phi == 0) + return a * (1-es); + + /* eq. 13b */ + if (fabs(phi) == M_PI_2) + return a / sqrt(1-es); + + /* eq. 13 */ + return (a * (1 - es) ) / pow(1 - es*s*s, 1.5); + +} + + +static LPZ calc_standard_params(LPZ lpz, PJ *P) { + struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; + double dphi, dlam, dh; + + /* sines and cosines */ + double slam = sin(lpz.lam); + double clam = cos(lpz.lam); + double sphi = sin(lpz.phi); + double cphi = cos(lpz.phi); + + /* ellipsoid parameters and differences */ + double f = P->f, a = P->a; + double dx = Q->dx, dy = Q->dy, dz = Q->dz; + double da = Q->da, df = Q->df; + + /* ellipsoid radii of curvature */ + double rho = RM(a, P->es, lpz.phi); + double nu = RN(a, P->e2s, lpz.phi); + + /* delta phi */ + dphi = (-dx*sphi*clam) - (dy*sphi*slam) + (dz*cphi) + + ((nu * P->es * sphi * cphi * da) / a) + + (sphi*cphi * ( rho/(1-f) + nu*(1-f))*df); + dphi /= (rho + lpz.z); + + /* delta lambda */ + dlam = (-dx*slam + dy*clam) / ((nu+lpz.z)*cphi); + + /* delta h */ + dh = dx*cphi*clam + dy*cphi*slam + dz*sphi - (a/nu)*da + nu*(1-f)*sphi*sphi*df; + + lpz.phi = dphi; + lpz.lam = dlam; + lpz.z = dh; + + return lpz; +} + + +static LPZ calc_abridged_params(LPZ lpz, PJ *P) { + struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; + double dphi, dlam, dh; + + /* sines and cosines */ + double slam = sin(lpz.lam); + double clam = cos(lpz.lam); + double sphi = sin(lpz.phi); + double cphi = cos(lpz.phi); + + /* ellipsoid parameters and differences */ + double dx = Q->dx, dy = Q->dy, dz = Q->dz; + double da = Q->da, df = Q->df; + double adffda = (P->a*df + P->f*da); + + /* delta phi */ + dphi = -dx*sphi*clam - dy*sphi*slam + dz*cphi + adffda*sin(2*lpz.phi); + dphi /= RM(P->a, P->es, lpz.phi); + + /* delta lambda */ + dlam = -dx*slam + dy*clam; + dlam /= RN(P->a, P->e2s, lpz.phi)*cphi; + + /* delta h */ + dh = dx*cphi*clam + dy*cphi*slam + dz*sphi - da + adffda*sphi*sphi; + + /* offset coordinate */ + lpz.phi = dphi; + lpz.lam = dlam; + lpz.z = dh; + + return lpz; +} + + +static XY forward_2d(LP lp, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + + point.lp = lp; + point.xyz = forward_3d(point.lpz, P); + + return point.xy; +} + + +static LP reverse_2d(XY xy, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + + point.xy = xy; + point.xyz.z = 0; + point.lpz = reverse_3d(point.xyz, P); + + return point.lp; +} + + +static XYZ forward_3d(LPZ lpz, PJ *P) { + struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + + point.lpz = lpz; + + /* calculate parameters depending on the mode we are in */ + if (Q->abridged) { + lpz = calc_abridged_params(lpz, P); + } else { + lpz = calc_standard_params(lpz, P); + } + + /* offset coordinate */ + point.lpz.phi += lpz.phi; + point.lpz.lam += lpz.lam; + point.lpz.z += lpz.z; + + return point.xyz; +} + + +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { + obs.xyz = forward_3d(obs.lpz, P); + return obs; +} + + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + LPZ lpz; + + /* calculate parameters depending on the mode we are in */ + point.xyz = xyz; + if (Q->abridged) + lpz = calc_abridged_params(point.lpz, P); + else + lpz = calc_standard_params(point.lpz, P); + + /* offset coordinate */ + point.lpz.phi -= lpz.phi; + point.lpz.lam -= lpz.lam; + point.lpz.z -= lpz.z; + + return point.lpz; +} + + +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { + obs.lpz = reverse_3d(obs.xyz, P); + return obs; +} + + +PJ *TRANSFORMATION(molodensky,1) { + int count_required_params = 0; + struct pj_opaque_molodensky *Q = static_cast(pj_calloc(1, sizeof(struct pj_opaque_molodensky))); + if (nullptr==Q) + return pj_default_destructor(P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + /* read args */ + if (pj_param(P->ctx, P->params, "tdx").i) { + count_required_params ++; + Q->dx = pj_param(P->ctx, P->params, "ddx").f; + } + + if (pj_param(P->ctx, P->params, "tdy").i) { + count_required_params ++; + Q->dy = pj_param(P->ctx, P->params, "ddy").f; + } + + if (pj_param(P->ctx, P->params, "tdz").i) { + count_required_params ++; + Q->dz = pj_param(P->ctx, P->params, "ddz").f; + } + + if (pj_param(P->ctx, P->params, "tda").i) { + count_required_params ++; + Q->da = pj_param(P->ctx, P->params, "dda").f; + } + + if (pj_param(P->ctx, P->params, "tdf").i) { + count_required_params ++; + Q->df = pj_param(P->ctx, P->params, "ddf").f; + } + + Q->abridged = pj_param(P->ctx, P->params, "tabridged").i; + + /* We want all parameters (except +abridged) to be set */ + if (count_required_params == 0) + return pj_default_destructor(P, PJD_ERR_NO_ARGS); + + if (count_required_params != 5) + return pj_default_destructor(P, PJD_ERR_MISSING_ARGS); + + return P; +} diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp new file mode 100644 index 00000000..b3da906d --- /dev/null +++ b/src/transformations/vgridshift.cpp @@ -0,0 +1,144 @@ +#define PJ_LIB__ + +#include +#include +#include +#include + +#include "proj_internal.h" +#include "projects.h" + +PROJ_HEAD(vgridshift, "Vertical grid shift"); + +namespace { // anonymous namespace +struct pj_opaque_vgridshift { + double t_final; + double t_epoch; + double forward_multiplier; +}; +} // anonymous namespace + +static XYZ forward_3d(LPZ lpz, PJ *P) { + struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + point.lpz = lpz; + + if (P->vgridlist_geoid != nullptr) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + point.xyz.z += Q->forward_multiplier * proj_vgrid_value(P, point.lp); + } + + return point.xyz; +} + + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + point.xyz = xyz; + + if (P->vgridlist_geoid != nullptr) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + point.xyz.z -= Q->forward_multiplier * proj_vgrid_value(P, point.lp); + } + + return point.lpz; +} + + +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { + struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; + PJ_COORD point = obs; + + /* If transformation is not time restricted, we always call it */ + if (Q->t_final==0 || Q->t_epoch==0) { + point.xyz = forward_3d (obs.lpz, P); + return point; + } + + /* Time restricted - only apply transform if within time bracket */ + if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) + point.xyz = forward_3d (obs.lpz, P); + + + return point; +} + +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { + struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque; + PJ_COORD point = obs; + + /* If transformation is not time restricted, we always call it */ + if (Q->t_final==0 || Q->t_epoch==0) { + point.lpz = reverse_3d (obs.xyz, P); + return point; + } + + /* Time restricted - only apply transform if within time bracket */ + if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) + point.lpz = reverse_3d (obs.xyz, P); + + return point; +} + + +PJ *TRANSFORMATION(vgridshift,0) { + struct pj_opaque_vgridshift *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_vgridshift))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = (void *) Q; + + if (!pj_param(P->ctx, P->params, "tgrids").i) { + proj_log_error(P, "vgridshift: +grids parameter missing."); + return pj_default_destructor(P, PJD_ERR_NO_ARGS); + } + + /* TODO: Refactor into shared function that can be used */ + /* by both vgridshift and hgridshift */ + if (pj_param(P->ctx, P->params, "tt_final").i) { + Q->t_final = pj_param (P->ctx, P->params, "dt_final").f; + if (Q->t_final == 0) { + /* a number wasn't passed to +t_final, let's see if it was "now" */ + /* and set the time accordingly. */ + if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) { + time_t now; + struct tm *date; + time(&now); + date = localtime(&now); + Q->t_final = 1900.0 + date->tm_year + date->tm_yday/365.0; + } + } + } + + if (pj_param(P->ctx, P->params, "tt_epoch").i) + Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f; + + /* historical: the forward direction subtracts the grid offset. */ + Q->forward_multiplier = -1.0; + if (pj_param(P->ctx, P->params, "tmultiplier").i) { + Q->forward_multiplier = pj_param(P->ctx, P->params, "dmultiplier").f; + } + + /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */ + proj_vgrid_init(P, "grids"); + + /* Was gridlist compiled properly? */ + if ( proj_errno(P) ) { + proj_log_error(P, "vgridshift: could not find required grid(s)."); + return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = nullptr; + P->inv = nullptr; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + return P; +} -- cgit v1.2.3