Merge pull request #588 from kbevers/kinematic-grids

Kinematic gridshifting

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