Merge pull request #588 from kbevers/kinematic-grids

Kinematic gridshifting

11 files changed, 734 insertions, 269 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/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 0adc9e00..659039ab 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -4,7 +4,6 @@
 
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
 
-
 static XYZ forward_3d(LPZ lpz, PJ *P) {
     PJ_TRIPLET point;
     point.lpz = lpz;
@@ -12,9 +11,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
     if (P->gridlist != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        pj_apply_gridshift_3( P->ctx, P->gridlist,
-                              P->gridlist_count, 1, 1, 0,
-                              &point.xyz.x, &point.xyz.y, &point.xyz.z );
+        point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD);
     }
 
     return point.xyz;
@@ -28,9 +25,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
     if (P->gridlist != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        pj_apply_gridshift_3( P->ctx, P->gridlist,
-                              P->gridlist_count, 0, 1, 0,
-                              &point.xyz.x, &point.xyz.y, &point.xyz.z );
+        point.lp = proj_hgrid_apply(P, point.lp, PJ_INV);
     }
 
     return point.lpz;
@@ -52,29 +47,8 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) {
 
 
 
-#if 0
-static XY forward_xy(LP lp, PJ *P) {
-    PJ_TRIPLET point;
-    point.lp = lp;
-    point.lpz.z = 0;
-    point.xyz = forward_3d (point.lpz, P);
-    return point.xy;
-}
-
-
-static LP reverse_lp(XY xy, PJ *P) {
-    PJ_TRIPLET point;
-    point.xy = xy;
-    point.xyz.z = 0;
-    point.lpz = reverse_3d (point.xyz, P);
-    return point.lp;
-}
-#endif
-
-
-
 PJ *PROJECTION(hgridshift) {
-    
+
     P->fwdobs = forward_obs;
     P->invobs = reverse_obs;
     P->fwd3d  = forward_3d;
@@ -90,14 +64,12 @@ PJ *PROJECTION(hgridshift) {
         return pj_default_destructor (P, PJD_ERR_NO_ARGS);
     }
 
-    /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */
-    P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s,
-                                             &(P->gridlist_count) );
 
+    proj_hgrid_init(P, "grids");
     /* Was gridlist compiled properly? */
-    if ( pj_ctx_get_errno(pj_get_ctx(P)) ) {
+    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 pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
     }
 
     return P;
@@ -119,26 +91,28 @@ int pj_hgridshift_selftest (void) {
         proj_destroy (P);
         return 99;
     }
-    
+
     /* fail on purpose: open non-existing grid */
     P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=@nonexistinggrid.gsb,anothernonexistinggrid.gsb");
     if (0!=P) {
         proj_destroy (P);
         return 999;
     }
-        
+
     /* Failure most likely means the grid is missing */
     P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80");
     if (0==P)
         return 10;
-        
+
     a = proj_obs_null;
     a.coo.lpz.lam = PJ_TORAD(173);
     a.coo.lpz.phi = PJ_TORAD(-45);
-    
+
     dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
-    if (dist > 0.00000001)
+    if (dist > 0.00000001) {
+        printf("dist: %f\n",dist);
         return 1;
+    }
 
     expect.coo.lpz.lam = PJ_TORAD(172.999892181021551);
     expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613);
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index ededd544..691f791b 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -9,14 +9,10 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
     PJ_TRIPLET point;
     point.lpz = lpz;
 
-    if (P->gridlist != NULL) {
+    if (P->vgridlist_geoid != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        pj_apply_vgridshift( P, "sgrids",
-                             &(P->gridlist),
-                             &(P->gridlist_count),
-                             1, 1, 0,
-                             &point.xyz.x, &point.xyz.y, &point.xyz.z );
+        point.xyz.z -= proj_vgrid_value(P, point.lp);
     }
 
     return point.xyz;
@@ -27,14 +23,10 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
     PJ_TRIPLET point;
     point.xyz = xyz;
 
-    if (P->gridlist != NULL) {
+    if (P->vgridlist_geoid != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        pj_apply_vgridshift( P, "sgrids",
-                             &(P->gridlist),
-                             &(P->gridlist_count),
-                             0, 1, 0,
-                             &point.xyz.x, &point.xyz.y, &point.xyz.z );
+        point.xyz.z += proj_vgrid_value(P, point.lp);
     }
 
     return point.lpz;
@@ -61,14 +53,13 @@ PJ *PROJECTION(vgridshift) {
         return pj_default_destructor(P, PJD_ERR_NO_ARGS);
     }
 
-    /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */
-    P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s,
-                                             &(P->gridlist_count) );
+    /* 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 ( pj_ctx_get_errno(P->ctx) ) {
+    if ( proj_errno(P) ) {
         proj_log_error(P, "vgridshift: could not find required grid(s).");
-        return pj_default_destructor(P, -38);
+        return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
     }
 
     P->fwdobs = forward_obs;
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..edb683b0 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -61,7 +61,9 @@ 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_apply_gridshift.c b/src/pj_apply_gridshift.c
index 91e2de26..45887abd 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -2,7 +2,7 @@
  * Project:  PROJ.4
  * Purpose:  Apply datum shifts based on grid shift files (normally NAD27 to
  *           NAD83 or the reverse).  This module is responsible for keeping
- *           a list of loaded grids, and calling with each one that is 
+ *           a list of loaded grids, and calling with each one that is
  *           allowed for a given datum (expressed as the nadgrids= parameter).
  * Author:   Frank Warmerdam, warmerdam@pobox.com
  *
@@ -30,9 +30,10 @@
 
 #define PJ_LIB__
 
-#include <projects.h>
 #include <string.h>
 #include <math.h>
+#include "proj_internal.h"
+#include "projects.h"
 
 /************************************************************************/
 /*                         pj_apply_gridshift()                         */
@@ -43,7 +44,7 @@
 /*      it to honour our public api.                                    */
 /************************************************************************/
 
-int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, 
+int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
                         long point_count, int point_offset,
                         double *x, double *y, double *z )
 
@@ -51,16 +52,16 @@ int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
     PJ_GRIDINFO **gridlist;
     int           grid_count;
     int           ret;
-    
+
     gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count );
 
     if( gridlist == NULL || grid_count == 0 )
         return ctx->last_errno;
 
-    ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse, 
+    ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse,
                                 point_count, point_offset, x, y, z );
 
-    /* 
+    /*
     ** Note this frees the array of grid list pointers, but not the grids
     ** which is as intended.  The grids themselves live on.
     */
@@ -72,13 +73,13 @@ int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
 /************************************************************************/
 /*                        pj_apply_gridshift_2()                        */
 /*                                                                      */
-/*      This implementation takes uses the gridlist from a coordinate    */
+/*      This implementation uses the gridlist from a coordinate         */
 /*      system definition.  If the gridlist has not yet been            */
 /*      populated in the coordinate system definition we set it up      */
 /*      now.                                                            */
 /************************************************************************/
 
-int pj_apply_gridshift_2( PJ *defn, int inverse, 
+int pj_apply_gridshift_2( PJ *defn, int inverse,
                           long point_count, int point_offset,
                           double *x, double *y, double *z )
 
@@ -86,10 +87,10 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
     if( defn->catalog_name != NULL )
         return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
                                       x, y, z );
-                                      
+
     if( defn->gridlist == NULL )
     {
-        defn->gridlist = 
+        defn->gridlist =
             pj_gridlist_from_nadgrids( pj_get_ctx( defn ),
                                        pj_param(defn->ctx, defn->params,"snadgrids").s,
                                        &(defn->gridlist_count) );
@@ -97,12 +98,76 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
         if( defn->gridlist == NULL || defn->gridlist_count == 0 )
             return defn->ctx->last_errno;
     }
-     
+
     return pj_apply_gridshift_3( pj_get_ctx( defn ),
-                                 defn->gridlist, defn->gridlist_count, inverse, 
+                                 defn->gridlist, defn->gridlist_count, inverse,
                                  point_count, point_offset, x, y, z );
 }
 
+/************************************************************************/
+/*                             find_ctable()                            */
+/*                                                                      */
+/*    Determine which grid is the correct given an input coordinate.    */
+/************************************************************************/
+
+static struct CTABLE* find_ctable(projCtx ctx, LP input, int grid_count, PJ_GRIDINFO **tables) {
+    int itable;
+    double epsilon;
+    struct CTABLE *ct = NULL;
+
+    /* keep trying till we find a table that works */
+    for( itable = 0; itable < grid_count; itable++ )
+    {
+
+        PJ_GRIDINFO *gi = tables[itable];
+        ct = gi->ct;
+        epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0;
+        /* skip tables that don't match our point at all.  */
+        if ( ct->ll.phi - epsilon > input.phi
+                || ct->ll.lam - epsilon > input.lam
+                || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon < input.phi)
+                || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon < input.lam) ) {
+            continue;
+        }
+
+        /* If we have child nodes, check to see if any of them apply. */
+        while( gi->child )
+        {
+            PJ_GRIDINFO *child;
+
+            for( child = gi->child; child != NULL; child = child->next )
+            {
+                struct CTABLE *ct1 = child->ct;
+                epsilon = (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0;
+
+                if( ct1->ll.phi - epsilon > input.phi
+                    || ct1->ll.lam - epsilon > input.lam
+                    || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon < input.phi)
+                    || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon < input.lam) ) {
+                    continue;
+                }
+                break;
+            }
+
+            /* If we didn't find a child then nothing more to do */
+            if( child == NULL ) break;
+
+            /* Otherwise use the child, first checking it's children */
+            gi = child;
+            ct = child->ct;
+        }
+        /* load the grid shift info if we don't have it. */
+        if( ct->cvs == NULL) {
+            if (!pj_gridinfo_load( ctx, gi ) ) {
+                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+                return NULL;
+            }
+        }
+        /* if we get this far we have found a suitable grid */
+        break;
+    }
+    return ct;
+}
 
 /************************************************************************/
 /*                        pj_apply_gridshift_3()                        */
@@ -113,16 +178,16 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
 int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
                           int inverse, long point_count, int point_offset,
                           double *x, double *y, double *z )
-
 {
     int  i;
+    struct CTABLE *ct;
     static int debug_count = 0;
     (void) z;
 
     if( tables == NULL || grid_count == 0 )
     {
-        pj_ctx_set_errno( ctx, -38);
-        return -38;
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return PJD_ERR_FAILED_TO_LOAD_GRID;
     }
 
     ctx->last_errno = 0;
@@ -138,100 +203,39 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
         output.phi = HUGE_VAL;
         output.lam = HUGE_VAL;
 
-        /* keep trying till we find a table that works */
-        for( itable = 0; itable < grid_count; itable++ )
-        {
-            PJ_GRIDINFO *gi = tables[itable];
-            struct CTABLE *ct = gi->ct;
-            double epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0;
-
-            /* skip tables that don't match our point at all.  */
-            if( ct->ll.phi - epsilon > input.phi 
-                || ct->ll.lam - epsilon > input.lam
-                || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon 
-                    < input.phi)
-                || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon 
-                    < input.lam) )
-                continue;
-
-            /* If we have child nodes, check to see if any of them apply. */
-            while( gi->child )
-            {
-                PJ_GRIDINFO *child;
-
-                for( child = gi->child; child != NULL; child = child->next )
-                {
-                    struct CTABLE *ct1 = child->ct;
-                    epsilon = 
-                        (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0;
-
-                    if( ct1->ll.phi - epsilon > input.phi 
-                        || ct1->ll.lam - epsilon > input.lam
-                        || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon 
-                            < input.phi)
-                        || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon 
-                            < input.lam) )
-                        continue;
-
-                    break;
-                }
-
-                /* If we didn't find a child then nothing more to do */
-
-                if( child == NULL ) break;
-
-                /* Otherwise use the child, first checking it's children */
-
-                gi = child;
-                ct = child->ct;
-            }
+        ct = find_ctable(ctx, input, grid_count, tables);
+        output = nad_cvt( input, inverse, ct );
 
-            /* load the grid shift info if we don't have it. */
-            if( ct->cvs == NULL && !pj_gridinfo_load( ctx, gi ) )
-            {
-                pj_ctx_set_errno( ctx, -38 );
-                return -38;
-            }
-            
-            output = nad_cvt( input, inverse, ct );
-            if( output.lam != HUGE_VAL )
-            {
-                if( debug_count++ < 20 )
-                    pj_log( ctx, PJ_LOG_DEBUG_MINOR,
-                            "pj_apply_gridshift(): used %s", ct->id );
-                break;
-            }
-        }
+        if ( output.lam != HUGE_VAL && debug_count++ < 20 )
+            pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id );
 
-        if( output.lam == HUGE_VAL )
+        if ( output.lam == HUGE_VAL )
         {
             if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
             {
                 pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
                     "pj_apply_gridshift(): failed to find a grid shift table for\n"
                     "                      location (%.7fdW,%.7fdN)",
-                    x[io] * RAD_TO_DEG, 
+                    x[io] * RAD_TO_DEG,
                     y[io] * RAD_TO_DEG );
                 for( itable = 0; itable < grid_count; itable++ )
                 {
                     PJ_GRIDINFO *gi = tables[itable];
                     if( itable == 0 )
-                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                                "   tried: %s", gi->gridname );
+                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, "   tried: %s", gi->gridname );
                     else
-                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                                ",%s", gi->gridname );
+                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, ",%s", gi->gridname );
                 }
             }
 
-            /* 
-             * We don't actually have any machinery currently to set the 
-             * following macro, so this is mostly kept here to make it clear 
-             * how we ought to operate if we wanted to make it super clear 
+            /*
+             * We don't actually have any machinery currently to set the
+             * following macro, so this is mostly kept here to make it clear
+             * how we ought to operate if we wanted to make it super clear
              * that an error has occurred when points are outside our available
-             * datum shift areas.  But if this is on, we will find that "low 
-             * value" points on the fringes of some datasets will completely 
-             * fail causing lots of problems when it is more or less ok to 
+             * datum shift areas.  But if this is on, we will find that "low
+             * value" points on the fringes of some datasets will completely
+             * fail causing lots of problems when it is more or less ok to
              * just not apply a datum shift.  So rather than deal with
              * that we just fallback to no shift. (see also bug #45).
              */
@@ -252,3 +256,92 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
     return 0;
 }
 
+/**********************************************/
+int proj_hgrid_init(PJ* P, const char *grids) {
+/**********************************************
+
+  Initizalize and populate list of horizontal
+  grids.
+
+    Takes a PJ-object and the plus-parameter
+    name that is used in the proj-string to
+    specify the grids to load, e.g. "+grids".
+    The + should be left out here.
+
+    Returns the number of loaded grids.
+
+***********************************************/
+
+    /* prepend "s" to the "grids" string to allow usage with pj_param */
+    char *sgrids = (char *) pj_malloc( (strlen(grids)+1) *sizeof(char) );
+    sprintf(sgrids, "%s%s", "s", grids);
+
+    if (P->gridlist == NULL) {
+        P->gridlist = pj_gridlist_from_nadgrids(
+            P->ctx,
+            pj_param(P->ctx, P->params, sgrids).s,
+            &(P->gridlist_count)
+        );
+
+        if( P->gridlist == NULL || P->gridlist_count == 0 ) {
+            pj_dealloc(sgrids);
+            return 0;
+        }
+    }
+
+    if (P->gridlist_count == 0) {
+        proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    }
+
+    pj_dealloc(sgrids);
+    return P->gridlist_count;
+}
+
+/********************************************/
+/*           proj_hgrid_value()             */
+/*                                          */
+/*    Return coordinate offset in grid      */
+/********************************************/
+LP proj_hgrid_value(PJ *P, LP lp) {
+    struct CTABLE *ct;
+    LP out;
+
+    ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
+
+    /* normalize input to ll origin */
+    lp.lam -= ct->ll.lam;
+    lp.phi -= ct->ll.phi;
+    lp.lam = adjlon(lp.lam - M_PI) + M_PI;
+
+    out = nad_intr(lp, ct);
+
+    if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) {
+        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
+    }
+
+    return out;
+}
+
+LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction) {
+    struct CTABLE *ct;
+    int inverse;
+    LP out;
+
+    out.lam = HUGE_VAL; out.phi = HUGE_VAL;
+
+    ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
+
+    if (ct == NULL || ct->cvs == NULL) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+        return out;
+    }
+
+    inverse = direction == PJ_FWD ? 0 : 1;
+    out = nad_cvt(lp, inverse, ct);
+
+    if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
+        pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+
+    return out;
+
+}
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
index 35047a19..3c7cc210 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -28,23 +28,120 @@
 
 #define PJ_LIB__
 
-#include <projects.h>
 #include <string.h>
 #include <math.h>
+#include "proj_internal.h"
+#include "projects.h"
+
+static double pj_read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) {
+    int  itable = 0;
+    double value = HUGE_VAL;
+    double grid_x, grid_y;
+    int    grid_ix, grid_iy;
+    int    grid_ix2, grid_iy2;
+    float  *cvs;
+    /* do not deal with NaN coordinates */
+    if( input.phi != input.phi || input.lam != input.lam )
+        itable = *gridlist_count_p;
+
+    /* keep trying till we find a table that works */
+    for ( ; itable < *gridlist_count_p; itable++ )
+    {
+        PJ_GRIDINFO *gi = tables[itable];
+
+        ct = gi->ct;
+
+        /* skip tables that don't match our point at all.  */
+        if( ct->ll.phi > input.phi || ct->ll.lam > input.lam
+            || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi
+            || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam )
+            continue;
+
+        /* If we have child nodes, check to see if any of them apply. */
+        while( gi->child != NULL )
+        {
+            PJ_GRIDINFO *child;
+
+            for( child = gi->child; child != NULL; child = child->next )
+            {
+                struct CTABLE *ct1 = child->ct;
+
+                if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam
+                  || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi
+                  || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam)
+                    continue;
+
+                break;
+            }
+
+            /* we didn't find a more refined child node to use, so go with current grid */
+            if( child == NULL )
+            {
+                break;
+            }
+
+            /* Otherwise let's try for childrens children .. */
+            gi = child;
+            ct = child->ct;
+        }
+
+        /* load the grid shift info if we don't have it. */
+        if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) )
+        {
+            pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+            return PJD_ERR_FAILED_TO_LOAD_GRID;
+        }
+
+    }
+
+    /* Interpolation a location within the grid */
+    grid_x = (input.lam - ct->ll.lam) / ct->del.lam;
+    grid_y = (input.phi - ct->ll.phi) / ct->del.phi;
+    grid_ix = (int) floor(grid_x);
+    grid_iy = (int) floor(grid_y);
+    grid_x -= grid_ix;
+    grid_y -= grid_iy;
+
+    grid_ix2 = grid_ix + 1;
+    if( grid_ix2 >= ct->lim.lam )
+        grid_ix2 = ct->lim.lam - 1;
+    grid_iy2 = grid_iy + 1;
+    if( grid_iy2 >= ct->lim.phi )
+        grid_iy2 = ct->lim.phi - 1;
+
+    cvs = (float *) ct->cvs;
+    value = cvs[grid_ix + grid_iy * ct->lim.lam]
+        * (1.0-grid_x) * (1.0-grid_y)
+        + cvs[grid_ix2 + grid_iy * ct->lim.lam]
+        * (grid_x) * (1.0-grid_y)
+        + cvs[grid_ix + grid_iy2 * ct->lim.lam]
+        * (1.0-grid_x) * (grid_y)
+        + cvs[grid_ix2 + grid_iy2 * ct->lim.lam]
+        * (grid_x) * (grid_y);
+
+    /* nodata?  */
+    /* GTX official nodata value if  -88.88880f, but some grids also */
+    /* use other  big values for nodata (e.g naptrans2008.gtx has */
+    /* nodata values like -2147479936), so test them too */
+    if( value > 1000 || value < -1000 || value == -88.88880f )
+        value = HUGE_VAL;
+
+
+    return value;
+}
 
 /************************************************************************/
 /*                        pj_apply_vgridshift()                         */
 /*                                                                      */
-/*      This implementation takes uses the gridlist from a coordinate    */
+/*      This implementation takes uses the gridlist from a coordinate   */
 /*      system definition.  If the gridlist has not yet been            */
 /*      populated in the coordinate system definition we set it up      */
 /*      now.                                                            */
 /************************************************************************/
-
 int pj_apply_vgridshift( PJ *defn, const char *listname,
-                         PJ_GRIDINFO ***gridlist_p, 
+                         PJ_GRIDINFO ***gridlist_p,
                          int *gridlist_count_p,
-                         int inverse, 
+                         int inverse,
                          long point_count, int point_offset,
                          double *x, double *y, double *z )
 
@@ -52,22 +149,23 @@ int pj_apply_vgridshift( PJ *defn, const char *listname,
     int  i;
     static int debug_count = 0;
     PJ_GRIDINFO **tables;
+    struct CTABLE ct;
 
     if( *gridlist_p == NULL )
     {
-        *gridlist_p = 
-            pj_gridlist_from_nadgrids( pj_get_ctx(defn), 
+        *gridlist_p =
+            pj_gridlist_from_nadgrids( pj_get_ctx(defn),
                                        pj_param(defn->ctx,defn->params,listname).s,
                                        gridlist_count_p );
 
         if( *gridlist_p == NULL || *gridlist_count_p == 0 )
             return defn->ctx->last_errno;
     }
-     
+
     if( *gridlist_count_p == 0 )
     {
-        pj_ctx_set_errno( defn->ctx, -38);
-        return -38;
+        pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return PJD_ERR_FAILED_TO_LOAD_GRID;
     }
 
     tables = *gridlist_p;
@@ -75,127 +173,37 @@ int pj_apply_vgridshift( PJ *defn, const char *listname,
 
     for( i = 0; i < point_count; i++ )
     {
+        double value;
         long io = i * point_offset;
         LP   input;
-        int  itable = 0;
-        double value = HUGE_VAL;
 
         input.phi = y[io];
         input.lam = x[io];
 
-        /* do not deal with NaN coordinates */
-        if( input.phi != input.phi || input.lam != input.lam )
-            itable = *gridlist_count_p;
+        value = pj_read_vgrid_value(defn, input, gridlist_count_p, tables, &ct);
 
-        /* keep trying till we find a table that works */
-        for( ; itable < *gridlist_count_p; itable++ )
+        if( inverse )
+            z[io] -= value;
+        else
+            z[io] += value;
+        if( value != HUGE_VAL )
         {
-            PJ_GRIDINFO *gi = tables[itable];
-            struct CTABLE *ct = gi->ct;
-            double grid_x, grid_y;
-            int    grid_ix, grid_iy;
-            int    grid_ix2, grid_iy2;
-            float  *cvs;
-
-            /* skip tables that don't match our point at all.  */
-            if( ct->ll.phi > input.phi || ct->ll.lam > input.lam
-                || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi
-                || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam )
-                continue;
-
-            /* If we have child nodes, check to see if any of them apply. */
-            while( gi->child != NULL )
-            {
-                PJ_GRIDINFO *child;
-
-                for( child = gi->child; child != NULL; child = child->next )
-                {
-                    struct CTABLE *ct1 = child->ct;
-
-                    if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam
-                      || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi
-                      || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam)
-                        continue;
-
-                    break;
-                }
-
-                /* we didn't find a more refined child node to use, so go with current grid */
-                if( child == NULL )
-                {
-                    break;
-                }
-
-                /* Otherwise let's try for childrens children .. */
-                gi = child;
-                ct = child->ct;
-            }
-
-            /* load the grid shift info if we don't have it. */
-            if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) )
-            {
-                pj_ctx_set_errno( defn->ctx, -38 );
-                return -38;
-            }
-
-            /* Interpolation a location within the grid */
-            grid_x = (input.lam - ct->ll.lam) / ct->del.lam;
-            grid_y = (input.phi - ct->ll.phi) / ct->del.phi;
-            grid_ix = (int) floor(grid_x);
-            grid_iy = (int) floor(grid_y);
-            grid_x -= grid_ix;
-            grid_y -= grid_iy;
-
-            grid_ix2 = grid_ix + 1;
-            if( grid_ix2 >= ct->lim.lam )
-                grid_ix2 = ct->lim.lam - 1;
-            grid_iy2 = grid_iy + 1;
-            if( grid_iy2 >= ct->lim.phi )
-                grid_iy2 = ct->lim.phi - 1;
-
-            cvs = (float *) ct->cvs;
-            value = cvs[grid_ix + grid_iy * ct->lim.lam]
-                * (1.0-grid_x) * (1.0-grid_y)
-                + cvs[grid_ix2 + grid_iy * ct->lim.lam]
-                * (grid_x) * (1.0-grid_y)
-                + cvs[grid_ix + grid_iy2 * ct->lim.lam]
-                * (1.0-grid_x) * (grid_y)
-                + cvs[grid_ix2 + grid_iy2 * ct->lim.lam]
-                * (grid_x) * (grid_y);
-
-            /* nodata?  */
-            /* GTX official nodata value if  -88.88880f, but some grids also */
-            /* use other  big values for nodata (e.g naptrans2008.gtx has */
-            /* nodata values like -2147479936), so test them too */
-            if( value > 1000 || value < -1000 || value == -88.88880f )
-                value = HUGE_VAL;
-            else
-            {
-                if( inverse )
-                    z[io] -= value;
-                else
-                    z[io] += value;
-            }
-                
-            if( value != HUGE_VAL )
-            {
-                if( debug_count++ < 20 )
-                    pj_log( defn->ctx, PJ_LOG_DEBUG_MINOR, 
-                            "pj_apply_gridshift(): used %s",
-                            ct->id );
+            if( debug_count++ < 20 ) {
+                proj_log_trace(defn, "pj_apply_gridshift(): used %s", ct.id);
                 break;
             }
         }
 
         if( value == HUGE_VAL )
         {
+            int itable;
             char gridlist[3000];
 
-            pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
-                    "pj_apply_vgridshift(): failed to find a grid shift table for\n"
-                    "                       location (%.7fdW,%.7fdN)",
-                    x[io] * RAD_TO_DEG, 
-                    y[io] * RAD_TO_DEG );
+            proj_log_debug(defn,
+                "pj_apply_vgridshift(): failed to find a grid shift table for\n"
+                "                       location (%.7fdW,%.7fdN)",
+                x[io] * RAD_TO_DEG,
+                y[io] * RAD_TO_DEG );
 
             gridlist[0] = '\0';
             for( itable = 0; itable < *gridlist_count_p; itable++ )
@@ -212,10 +220,10 @@ int pj_apply_vgridshift( PJ *defn, const char *listname,
                 else
                     sprintf( gridlist+strlen(gridlist), ",%s", gi->gridname );
             }
-            pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
-                    "%s", gridlist );
-                
+
+            proj_log_debug(defn, "%s", gridlist);
             pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA );
+
             return PJD_ERR_GRID_AREA;
         }
     }
@@ -223,3 +231,64 @@ int pj_apply_vgridshift( PJ *defn, const char *listname,
     return 0;
 }
 
+/**********************************************/
+int proj_vgrid_init(PJ* P, const char *grids) {
+/**********************************************
+
+  Initizalize and populate gridlist.
+
+    Takes a PJ-object and the plus-parameter
+    name that is used in the proj-string to
+    specify the grids to load, e.g. "+grids".
+    The + should be left out here.
+
+    Returns the number of loaded grids.
+
+***********************************************/
+
+    /* prepend "s" to the "grids" string to allow usage with pj_param */
+    char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) );
+    sprintf(sgrids, "%s%s", "s", grids);
+
+    if (P->vgridlist_geoid == NULL) {
+        P->vgridlist_geoid = pj_gridlist_from_nadgrids(
+            P->ctx,
+            pj_param(P->ctx, P->params, sgrids).s,
+            &(P->vgridlist_geoid_count)
+        );
+
+        if( P->vgridlist_geoid == NULL || P->vgridlist_geoid_count == 0 ) {
+            pj_dealloc(sgrids);
+            return 0;
+        }
+    }
+
+    if (P->vgridlist_geoid_count == 0) {
+        proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    }
+
+    pj_dealloc(sgrids);
+    return P->vgridlist_geoid_count;
+}
+
+/***********************************************/
+double proj_vgrid_value(PJ *P, LP lp){
+/***********************************************
+
+  Read grid value at position lp in grids loaded
+  with proj_grid_init.
+
+  Returns the grid value of the given coordinate.
+
+************************************************/
+
+    struct CTABLE used_grid;
+    double value;
+    memset(&used_grid, 0, sizeof(struct CTABLE));
+
+    value = pj_read_vgrid_value(P, lp, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid);
+    proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value);
+
+    return value;
+}
+
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 5773637c..ba542dea 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -89,7 +89,12 @@ PJ_OBS pj_invobs (PJ_OBS obs, PJ *P);
 PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P);
 PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P);
 
-
+/* Grid functionality */
+int             proj_vgrid_init(PJ *P, const char *grids);
+int             proj_hgrid_init(PJ *P, const char *grids);
+double          proj_vgrid_value(PJ *P, LP lp);
+LP              proj_hgrid_value(PJ *P, LP lp);
+LP              proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction);
 
 /* High level functionality for handling thread contexts */
 enum proj_log_level {
@@ -112,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);