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.

1 files changed, 184 insertions, 91 deletions

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;
+
+}