Addition of 'deformation': Kinematic grid shifting.

Kinematic deformation models are used in some geodetic transformations.
This commit introduces the ability to do transformations involving a
gridded deformation/velocity model. For practical reasons a gridded deformation
model needs to be split into two seperate files, one for the horizontal
components and one for the vertical component. For this we use formats
already known to PROJ.4, namely the CTable/CTable2 and the GTX formats.
Grids are specified in the proj-string with +xy_grids and +z_grids.
Grid values are expected to be in m/year. The kinematic part of the
operation is controlled by the +t_epoch parameter, which is the central
epoch of the transformation. An observation epoch is also needed. It can
be specified either in the PJ_OBS input as the fourth element in the
coordinate, or in the proj-string with +t_obs. If +t_obs is present in
the proj-string it takes presedence over the value in the PJ_OBS
coordinate.

9 files changed, 410 insertions, 4 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 363f2cce..0bfcfb24 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -87,7 +87,7 @@ libproj_la_SOURCES = \
 	\
 	proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \
 	PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \
-	pj_internal.c
+	PJ_deformation.c pj_internal.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c
new file mode 100644
index 00000000..05858900
--- /dev/null
+++ b/src/PJ_deformation.c
@@ -0,0 +1,291 @@
+/***********************************************************************
+
+        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
+
+
+Corrections are done in cartesian space. Coordinates of the
+gridded model are in lat/long space and the corrections in the model are
+in cartesian space. Hence the input coordinates needs to be converted
+to lat/long space when searching for corrections in the grid. The
+corrections are then applied to the input coordinates 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. Both grids are expected to contain grid-values in units of
+m/year.
+
+************************************************************************
+* 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 <errno.h>
+#include <proj.h>
+#include "proj_internal.h"
+#include "projects.h"
+
+PROJ_HEAD(deformation, "Kinematic grid shift");
+
+#define TOL 1e-8
+#define MAX_ITERATIONS 10
+
+struct pj_opaque {
+    double t_obs;
+    double t_epoch;
+    int has_xy_grids, has_z_grids;
+    PJ *cart;
+};
+
+static XYZ get_xygrid_shift(PJ* P, XYZ cartesian) {
+    PJ_COORD geodetic, shift;
+
+    geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart);
+    shift.lp = proj_hgrid_value(P, geodetic.lp);
+
+    return shift.xyz;
+}
+
+static double get_zgrid_shift(PJ* P, XYZ cartesian) {
+    PJ_COORD geodetic;
+
+    geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart);
+
+    return proj_vgrid_value(P, geodetic.lp);
+}
+
+static XYZ reverse_hshift(PJ *P, XYZ input, double dt) {
+    XYZ out, delta, dif;
+    int i = MAX_ITERATIONS;
+
+    delta = get_xygrid_shift(P, input);
+
+    out.x = input.x + dt*delta.x;
+    out.y = input.y - dt*delta.y;
+    out.z = input.z;
+
+    do {
+        delta = get_xygrid_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;
+        out.x += dif.x;
+        out.y += dif.y;
+
+    } while ( --i && hypot(dif.x, dif.y) > TOL );
+
+    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_obs - Q->t_epoch;
+    } 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;
+    }
+
+    if (Q->has_xy_grids) {
+        shift = get_xygrid_shift(P, in.xyz);
+
+        out.xyz.x += dt * shift.x;
+        out.xyz.y += dt * shift.y;
+    }
+
+    if (Q->has_z_grids)
+        out.xyz.z += dt * get_zgrid_shift(P, in.xyz);
+
+    return out.xyz;
+}
+
+
+static PJ_OBS forward_obs(PJ_OBS in, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    double dt;
+    XYZ shift;
+    PJ_OBS out = in;
+
+    if (Q->t_obs != HUGE_VAL) {
+            dt = Q->t_obs - Q->t_epoch;
+        } else {
+            dt = in.coo.xyzt.t - Q->t_epoch;
+    }
+
+    if (Q->has_xy_grids) {
+        shift = get_xygrid_shift(P, in.coo.xyz);
+        out.coo.xyz.x += dt * shift.x;
+        out.coo.xyz.y += dt * shift.y;
+    }
+
+    if (Q->has_z_grids)
+        out.coo.xyz.z += dt * get_zgrid_shift(P, in.coo.xyz);
+
+    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;
+    }
+
+    if (Q->has_xy_grids) {
+        out.xyz = reverse_hshift(P, in, dt);
+    }
+
+    if (Q->has_z_grids)
+        out.xyz.z = in.z + dt * get_zgrid_shift(P, in);
+
+    return out.lpz;
+}
+
+static PJ_OBS reverse_obs(PJ_OBS in, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    PJ_OBS out = in;
+    double dt;
+
+    if (Q->t_obs != HUGE_VAL) {
+            dt = Q->t_epoch - Q->t_obs;
+        } else {
+            dt = Q->t_epoch - in.coo.xyzt.t;
+    }
+
+    if (Q->has_xy_grids)
+        out.coo.xyz = reverse_hshift(P, in.coo.xyz, dt);
+
+    if (Q->has_z_grids)
+        out.coo.xyz.z = in.coo.xyz.z + dt * get_zgrid_shift(P, in.coo.xyz);
+
+    return out;
+}
+
+static void *destructor(PJ *P, int errlev) {
+    if (0==P)
+        return 0;
+
+    if (0==P->opaque)
+        return pj_default_destructor (P, errlev);
+
+    if (P->opaque->cart)
+        P->opaque->cart->destructor (P->opaque->cart, errlev);
+
+    return pj_default_destructor(P, errlev);
+}
+
+
+PJ *PROJECTION(deformation) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return pj_default_destructor(P, ENOMEM);
+    P->opaque = (void *) Q;
+
+    Q->cart = proj_create(P->ctx, "+proj=cart");
+    if (Q->cart == 0)
+        return destructor(P, ENOMEM);
+
+    /* inherit ellipsoid definition from P to Q->cart (simpler than guessing */
+    /* how the ellipsoid was specified in original definition)               */
+    pj_inherit_ellipsoid_defs(P, Q->cart);
+
+    Q->has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i;
+    Q->has_z_grids  = pj_param(P->ctx, P->params, "tz_grids").i;
+
+    /* Build gridlists. P->gridlist and P->vgridlist_geoid can be empty if */
+    /* +xy_grids or +z_grids only ask for optional grids                   */
+    if (!Q->has_xy_grids && !Q->has_z_grids) {
+        proj_log_error(P, "deformation: At least one of either +xy_grids or +z_grids should be specified.");
+        return destructor(P, PJD_ERR_NO_ARGS );
+    }
+
+    if (Q->has_xy_grids) {
+        Q->has_xy_grids = proj_hgrid_init(P, "xy_grids");
+        if (proj_errno(P)) {
+            proj_log_error(P, "deformation: could not find requested grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
+    }
+
+    if (Q->has_z_grids) {
+        Q->has_z_grids = proj_vgrid_init(P, "z_grids");
+        if (proj_errno(P)) {
+            proj_log_error(P, "deformation: could not find requested 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->fwdobs = forward_obs;
+    P->invobs = reverse_obs;
+    P->fwd3d  = forward_3d;
+    P->inv3d  = reverse_3d;
+    P->fwd    = 0;
+    P->inv    = 0;
+
+    P->left  = PJ_IO_UNITS_METERS;
+    P->right = PJ_IO_UNITS_METERS;
+    P->destructor = destructor;
+
+    return P;
+}
+
+int pj_deformation_selftest (void) {return 10000;}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 053e9ef6..969ad64b 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -58,6 +58,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_collg.c
         PJ_comill.c
         PJ_crast.c
+        PJ_deformation.c
         PJ_denoy.c
         PJ_eck1.c
         PJ_eck2.c
diff --git a/src/makefile.vc b/src/makefile.vc
index 1330e9bb..15fd6428 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -61,7 +61,8 @@ support = \
 
 pipeline = \
 	proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
-	PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj
+	PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \
+	PJ_deformation.obj
 
 geodesic = geodesic.obj
 
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 72892ac2..d43feecb 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,12 +1,49 @@
 /* set ellipsoid parameters a and es */
-#include <projects.h>
 #include <string.h>
+#include "proj_internal.h"
+#include "projects.h"
 #define SIXTH .1666666666666666667  /* 1/6 */
 #define RA4   .04722222222222222222 /* 17/360 */
 #define RA6   .02215608465608465608 /* 67/3024 */
 #define RV4   .06944444444444444444 /* 5/72 */
 #define RV6   .04243827160493827160 /* 55/1296 */
 
+/* copy ellipsoidal parameters from src to dst */
+void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
+
+    /* The linear parameters */
+    dst->a  = src->a;
+    dst->b  = src->b;
+    dst->ra = src->ra;
+    dst->rb = src->rb;
+
+    /* The eccentricities */
+    dst->alpha   = src->alpha;
+    dst->e       = src->e;
+    dst->es      = src->es;
+    dst->e2      = src->e2;
+    dst->e2s     = src->e2s;
+    dst->e3      = src->e3;
+    dst->e3s     = src->e3s;
+    dst->one_es  = src->one_es;
+    dst->rone_es = src->rone_es;
+
+    /* The flattenings */
+    dst->f   = src->f;
+    dst->f2  = src->f2;
+    dst->n   = src->n;
+    dst->rf  = src->rf;
+    dst->rf2 = src->rf2;
+    dst->rn  = src->rn;
+
+    /* This one's for GRS80 */
+    dst->J = src->J;
+
+    /* es and a before any +proj related adjustment */
+    dst->es_orig = src->es_orig;
+    dst->a_orig  = src->a_orig;
+}
+
 /* initialize geographic shape parameters */
 int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
     int i;
diff --git a/src/pj_list.h b/src/pj_list.h
index bf287219..76c37126 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -25,6 +25,7 @@ PROJ_HEAD(chamb, "Chamberlin Trimetric")
 PROJ_HEAD(collg, "Collignon")
 PROJ_HEAD(comill, "Compact Miller")
 PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)")
+PROJ_HEAD(deformation, "Kinematic grid shift")
 PROJ_HEAD(denoy, "Denoyer Semi-Elliptical")
 PROJ_HEAD(eck1, "Eckert I")
 PROJ_HEAD(eck2, "Eckert II")
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 29e73e70..ba542dea 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -117,6 +117,7 @@ void proj_log_trace (PJ *P, const char *fmt, ...);
 /*void proj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));*/
 void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log);
 
+void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst);
 
 /* Lowest level: Minimum support for fileapi */
 void proj_fileapi_set (PJ *P, void *fileapi);
diff --git a/test/gie/deformation.gie b/test/gie/deformation.gie
new file mode 100644
index 00000000..74a6b25d
--- /dev/null
+++ b/test/gie/deformation.gie
@@ -0,0 +1,73 @@
+BEGIN
+===============================================================================
+Test for the deformation operation - Kinematic Gridshifting
+
+For all the deformation tests the alaska and egm96_15.gtx grids are used even
+though they are not parts of a deformation model, they are in the proper format
+and for testing purposes it doesn't really matter all that much...
+
+The input coordinate is located at lon=60, lam=-160 - somewhere in Alaska.
+
+===============================================================================
+
+
+-------------------------------------------------------------------------------
+Test using only horizontal grid and +tobs parameter
+-------------------------------------------------------------------------------
+OPERATION   +proj=deformation +xy_grids=alaska +t_epoch=2016.0 +t_obs=2000.0
+-------------------------------------------------------------------------------
+TOLERANCE   0.000001 m
+ACCEPT      -3004295.5882503074 -1093474.1690603832 5500477.1338251457
+EXPECT      -3004295.5888766116 -1093474.1688513425 5500477.1338251457
+ROUNDTRIP   1000
+
+-------------------------------------------------------------------------------
+Test using only vertical grid and +tobs parameter
+-------------------------------------------------------------------------------
+OPERATION   +proj=deformation +z_grids=egm96_15.gtx +t_epoch=2016.0 +t_obs=2000.0
+-------------------------------------------------------------------------------
+TOLERANCE   0.000001 m
+ACCEPT      -3004295.5882503074 -1093474.1690603832 5500477.1338251457
+EXPECT      -3004295.5882503074 -1093474.1690603832 5500234.008855661
+ROUNDTRIP   1000
+
+-------------------------------------------------------------------------------
+Test using both horizontal and vertical grids as well as the +tobs parameter
+-------------------------------------------------------------------------------
+OPERATION   +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=2016.0 +t_obs=2000.0
+-------------------------------------------------------------------------------
+TOLERANCE   0.000001 m
+ACCEPT      -3004295.5882503074 -1093474.1690603832 5500477.1338251457
+EXPECT      -3004295.5888766116 -1093474.1688513425 5500234.008855661
+ROUNDTRIP   1000
+
+-------------------------------------------------------------------------------
+Test using only horizontal grid
+-------------------------------------------------------------------------------
+OPERATION   +proj=deformation +xy_grids=alaska +t_epoch=2016.0
+-------------------------------------------------------------------------------
+TOLERANCE   0.000001 m
+ACCEPT      -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
+EXPECT      -3004295.5888766116 -1093474.1688513425 5500477.1338251457 2000.0
+ROUNDTRIP   1000
+
+-------------------------------------------------------------------------------
+Test using only vertical grid
+-------------------------------------------------------------------------------
+OPERATION   +proj=deformation +z_grids=egm96_15.gtx +t_epoch=2016.0
+-------------------------------------------------------------------------------
+TOLERANCE   0.000001 m
+ACCEPT      -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
+EXPECT      -3004295.5882503074 -1093474.1690603832 5500234.008855661  2000.0
+ROUNDTRIP   1000
+
+-------------------------------------------------------------------------------
+Test using both horizontal and vertical grids
+-------------------------------------------------------------------------------
+OPERATION   +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=2016.0
+-------------------------------------------------------------------------------
+TOLERANCE   0.000001 m
+ACCEPT      -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
+EXPECT      -3004295.5888766116 -1093474.1688513425 5500234.008855661  2000.0
+ROUNDTRIP   1000
+END
diff --git a/travis/install.sh b/travis/install.sh
index 1247435b..62c2677f 100755
--- a/travis/install.sh
+++ b/travis/install.sh
@@ -81,7 +81,8 @@ if [ $TRAVIS_OS_NAME == "osx" ]; then
 make -j3
 make check
 PROJ_LIB=$GRIDDIR ./src/proj -VC
-./src/gie ./test/gie/builtins.gie
+PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/builtins.gie
+PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/deformation.gie
 
 # install & run the working GIGS test
   # create locations that pyproj understands