Refactor grid shift functions.

This refactoring of the grid shift functions allows for easier access to
the actual grid values, as well as making it possible to perform
horizontal grid shift on a single coordinate without using the
pj_apply_gridshift* functions. The latter simplifies the execution path
of the forward and inverse functions in PJ_hgridshift.c.

This commit introduces proj_*grid_init, proj_*grid_value
and proj_hgrid_apply. The init functions initialises horizontal and
vertical grids respectivelive (wrappers for pj_gridlist_from_nadgrids
with simpler parameters). The proj_*grid_value functions returns the
specific grid value at coordinate lp. The proj_hgrid_apply function
applies the grid offset to the input coordinate and outputs the adjusted
coordinate.

5 files changed, 398 insertions, 266 deletions

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/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/proj_internal.h b/src/proj_internal.h
index 5773637c..29e73e70 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 {