horizontal grid shift: rework to no longer load whole grid into memory

1 files changed, 260 insertions, 242 deletions

diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
index 2f46ca11..5a5404da 100644
--- a/src/apply_gridshift.cpp
+++ b/src/apply_gridshift.cpp
@@ -28,7 +28,9 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#define PJ_LIB__
+#ifndef FROM_PROJ_CPP
+#define FROM_PROJ_CPP
+#endif
 
 #include <math.h>
 #include <stdio.h>
@@ -36,231 +38,62 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj/internal/internal.hpp"
 
-/************************************************************************/
-/*                         pj_apply_gridshift()                         */
-/*                                                                      */
-/*      This is the externally callable interface - part of the         */
-/*      public API - though it is not used internally any more and I    */
-/*      doubt it is used by any other applications.  But we preserve    */
-/*      it to honour our public api.                                    */
-/************************************************************************/
+using namespace NS_PROJ;
 
-int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
-                        long point_count, int point_offset,
-                        double *x, double *y, double *z )
+// ---------------------------------------------------------------------------
 
+static const HorizontalShiftGrid* findGrid(
+    const std::vector<std::unique_ptr<HorizontalShiftGridSet>>& grids,
+    PJ_LP input)
 {
-    PJ_GRIDINFO **gridlist;
-    int           grid_count;
-    int           ret;
-
-    gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count );
-
-    if( gridlist == nullptr || grid_count == 0 )
-    {
-        pj_dalloc( gridlist );
-        return ctx->last_errno;
-    }
-
-    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.
-    */
-    pj_dalloc( gridlist );
-
-    return ret;
-}
-
-/************************************************************************/
-/*                        pj_apply_gridshift_2()                        */
-/*                                                                      */
-/*      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,
-                          long point_count, int point_offset,
-                          double *x, double *y, double *z )
-
-{
-    if( defn->gridlist == nullptr )
-    {
-        defn->gridlist =
-            pj_gridlist_from_nadgrids( pj_get_ctx( defn ),
-                                       pj_param(defn->ctx, defn->params,"snadgrids").s,
-                                       &(defn->gridlist_count) );
-
-        if( defn->gridlist == nullptr || defn->gridlist_count == 0 )
-            return defn->ctx->last_errno;
-    }
-
-    return pj_apply_gridshift_3( pj_get_ctx( defn ),
-                                 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, PJ_LP input, int grid_count, PJ_GRIDINFO **tables) {
-    int itable;
-
-    /* keep trying till we find a table that works */
-    for( itable = 0; itable < grid_count; itable++ )
+    for( const auto& gridset: grids )
     {
-
-        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 != nullptr; 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 == nullptr ) 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 == nullptr) {
-            if (!pj_gridinfo_load( ctx, gi ) ) {
-                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                return nullptr;
-            }
-        }
-        /* if we get this far we have found a suitable grid */
-        return ct;
+        auto grid = gridset->gridAt(input.lam, input.phi);
+        if( grid )
+            return grid;
     }
-
     return nullptr;
 }
 
-/************************************************************************/
-/*                        pj_apply_gridshift_3()                        */
-/*                                                                      */
-/*      This is the real workhorse, given a gridlist.                   */
-/************************************************************************/
+// ---------------------------------------------------------------------------
 
-int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_count,
-                          int inverse, long point_count, int point_offset,
-                          double *x, double *y, double *z )
+static std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>>
+                                        getListOfGridSets(PJ_CONTEXT* ctx,
+                                                          const char* grids)
 {
-    int  i;
-    struct CTABLE *ct;
-    static int debug_count = 0;
-    (void) z;
-
-    if( gridlist== nullptr || gridlist_count == 0 )
+    std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>> list;
+    auto listOfGrids = internal::split(std::string(grids), ',');
+    for( const auto& grid: listOfGrids )
     {
-        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
-        return PJD_ERR_FAILED_TO_LOAD_GRID;
-    }
-
-    ctx->last_errno = 0;
-
-    for( i = 0; i < point_count; i++ )
-    {
-        long io = i * point_offset;
-        PJ_LP   input, output;
-        int  itable;
-
-        input.phi = y[io];
-        input.lam = x[io];
-        output.phi = HUGE_VAL;
-        output.lam = HUGE_VAL;
-
-        ct = find_ctable(ctx, input, gridlist_count, gridlist);
-        if( ct != nullptr )
+        const char* gridname = grid.c_str();
+        bool canFail = false;
+        if( gridname[0] == '@' )
         {
-            output = nad_cvt( input, inverse, ct );
-
-            if ( output.lam != HUGE_VAL && debug_count++ < 20 )
-                pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id );
+            canFail = true;
+            gridname ++;
         }
-
-        if ( output.lam == HUGE_VAL )
+        auto gridSet = HorizontalShiftGridSet::open(ctx, gridname);
+        if( !gridSet )
         {
-            if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
+            if( !canFail )
             {
-                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,
-                    y[io] * RAD_TO_DEG );
-                for( itable = 0; itable < gridlist_count; itable++ )
-                {
-                    PJ_GRIDINFO *gi = gridlist[itable];
-                    if( itable == 0 )
-                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, "   tried: %s", gi->gridname );
-                    else
-                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, ",%s", gi->gridname );
-                }
+                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+                return {};
             }
-
-            /*
-             * 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
-             * just not apply a datum shift.  So rather than deal with
-             * that we just fallback to no shift. (see also bug #45).
-             */
-#ifdef ERR_GRID_AREA_TRANSIENT_SEVERE
-            y[io] = HUGE_VAL;
-            x[io] = HUGE_VAL;
-#else
-            /* leave x/y unshifted. */
-#endif
         }
         else
         {
-            y[io] = output.phi;
-            x[io] = output.lam;
+            list.emplace_back(std::move(gridSet));
         }
     }
-
-    return 0;
+    return list;
 }
 
+
 /**********************************************/
-int proj_hgrid_init(PJ* P, const char *grids) {
+int proj_hgrid_init(PJ* P, const char *gridkey) {
 /**********************************************
 
   Initizalize and populate list of horizontal
@@ -275,29 +108,158 @@ int proj_hgrid_init(PJ* P, const char *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);
+    std::string key("s");
+    key += gridkey;
+    const char* grids = pj_param(P->ctx, P->params, key.c_str()).s;
+    if( grids == nullptr )
+        return 0;
 
-    if (P->gridlist == nullptr) {
-        P->gridlist = pj_gridlist_from_nadgrids(
-            P->ctx,
-            pj_param(P->ctx, P->params, sgrids).s,
-            &(P->gridlist_count)
-        );
+    P->hgrids = getListOfGridSets(P->ctx, grids);
+    return static_cast<int>(P->hgrids.size());
+}
 
-        if( P->gridlist == nullptr || P->gridlist_count == 0 ) {
-            pj_dealloc(sgrids);
-            return 0;
+typedef struct { pj_int32 lam, phi; } ILP;
+
+static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid) {
+	PJ_LP val, frct;
+	ILP indx;
+	int in;
+
+        const auto& extent = grid->extentAndRes();
+	t.lam /= extent.resLon;
+	indx.lam = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam));
+	t.phi /= extent.resLat;
+	indx.phi = isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi));
+
+	frct.lam = t.lam - indx.lam;
+	frct.phi = t.phi - indx.phi;
+	val.lam = val.phi = HUGE_VAL;
+	if (indx.lam < 0) {
+		if (indx.lam == -1 && frct.lam > 0.99999999999) {
+			++indx.lam;
+			frct.lam = 0.;
+		} else
+			return val;
+	} else if ((in = indx.lam + 1) >= grid->width()) {
+		if (in == grid->width() && frct.lam < 1e-11) {
+			--indx.lam;
+			frct.lam = 1.;
+		} else
+			return val;
+	}
+	if (indx.phi < 0) {
+		if (indx.phi == -1 && frct.phi > 0.99999999999) {
+			++indx.phi;
+			frct.phi = 0.;
+		} else
+			return val;
+	} else if ((in = indx.phi + 1) >= grid->height()) {
+		if (in == grid->height() && frct.phi < 1e-11) {
+			--indx.phi;
+			frct.phi = 1.;
+		} else
+			return val;
+	}
+	
+	float f00Lon = 0, f00Lat = 0;
+        float f10Lon = 0, f10Lat = 0;
+        float f01Lon = 0, f01Lat = 0;
+        float f11Lon = 0, f11Lat = 0;
+	if( !grid->valueAt(indx.lam, indx.phi, f00Lon, f00Lat) ||
+            !grid->valueAt(indx.lam + 1, indx.phi, f10Lon, f10Lat) ||
+            !grid->valueAt(indx.lam, indx.phi + 1, f01Lon, f01Lat) ||
+            !grid->valueAt(indx.lam + 1, indx.phi + 1, f11Lon, f11Lat) )
+        {
+            return val;
         }
+
+        double m10 = frct.lam;
+	double m11 = m10;
+        double m01 = 1. - frct.lam;
+	double m00 = m01;
+	m11 *= frct.phi;
+	m01 *= frct.phi;
+	frct.phi = 1. - frct.phi;
+	m00 *= frct.phi;
+	m10 *= frct.phi;
+	val.lam = m00 * f00Lon + m10 * f10Lon +
+			  m01 * f01Lon + m11 * f11Lon;
+	val.phi = m00 * f00Lat + m10 * f10Lat +
+			  m01 * f01Lat + m11 * f11Lat;
+	return val;
+}
+
+
+#define MAX_ITERATIONS 10
+#define TOL 1e-12
+
+static
+PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) {
+    PJ_LP t, tb,del, dif;
+    int i = MAX_ITERATIONS;
+    const double toltol = TOL*TOL;
+
+    if (in.lam == HUGE_VAL)
+        return in;
+
+    /* normalize input to ll origin */
+    tb = in;
+    const auto& extent = grid->extentAndRes();
+    tb.lam -= extent.westLon;
+    tb.phi -= extent.southLat;
+
+    tb.lam = adjlon (tb.lam - M_PI) + M_PI;
+
+    t = nad_intr (tb, grid);
+    if (t.lam == HUGE_VAL)
+        return t;
+
+    if (!inverse) {
+        in.lam -= t.lam;
+        in.phi += t.phi;
+        return in;
     }
 
-    if (P->gridlist_count == 0) {
-        proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    t.lam = tb.lam + t.lam;
+    t.phi = tb.phi - t.phi;
+
+    do {
+        del = nad_intr(t, grid);
+
+        /* This case used to return failure, but I have
+           changed it to return the first order approximation
+           of the inverse shift.  This avoids cases where the
+           grid shift *into* this grid came from another grid.
+           While we aren't returning optimally correct results
+           I feel a close result in this case is better than
+           no result.  NFW
+           To demonstrate use -112.5839956 49.4914451 against
+           the NTv2 grid shift file from Canada. */
+        if (del.lam == HUGE_VAL)
+            break;
+
+        dif.lam = t.lam - del.lam - tb.lam;
+           dif.phi = t.phi + del.phi - tb.phi;
+        t.lam -= dif.lam;
+        t.phi -= dif.phi;
+
+    } while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); /* prob. slightly faster than hypot() */
+
+    if (i==0) {
+        /* If we had access to a context, this should go through pj_log, and we should set ctx->errno */
+        if (getenv ("PROJ_DEBUG"))
+            fprintf( stderr, "Inverse grid shift iterator failed to converge.\n" );
+        t.lam = t.phi = HUGE_VAL;
+        return t;
     }
 
-    pj_dealloc(sgrids);
-    return P->gridlist_count;
+    /* and again: pj_log and ctx->errno */
+    if (del.lam==HUGE_VAL && getenv ("PROJ_DEBUG"))
+        fprintf (stderr, "Inverse grid shift iteration failed, presumably at grid edge.\nUsing first approximation.\n");
+
+    in.lam = adjlon (t.lam + extent.westLon);
+    in.phi = t.phi + extent.southLat;
+    return in;
 }
 
 /********************************************/
@@ -306,22 +268,22 @@ int proj_hgrid_init(PJ* P, const char *grids) {
 /*    Return coordinate offset in grid      */
 /********************************************/
 PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) {
-    struct CTABLE *ct;
     PJ_LP out = proj_coord_error().lp;
 
-    ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
-    if (ct == nullptr) {
-        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA);
+    const auto grid = findGrid(P->hgrids, lp);
+    if( !grid ) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA );
         return out;
     }
 
     /* normalize input to ll origin */
-    lp.lam -= ct->ll.lam;
-    lp.phi -= ct->ll.phi;
+    const auto& extent = grid->extentAndRes();
+    lp.lam -= extent.westLon;
+    lp.phi -= extent.southLat;
 
     lp.lam = adjlon(lp.lam - M_PI) + M_PI;
 
-    out = nad_intr(lp, ct);
+    out = nad_intr(lp, grid);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) {
         pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
@@ -330,33 +292,89 @@ PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) {
     return out;
 }
 
-PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) {
-    struct CTABLE *ct;
-    int inverse;
+static
+PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx,
+                       PJ_LP lp,
+                       PJ_DIRECTION direction,
+                       const std::vector<std::unique_ptr<HorizontalShiftGridSet>>& grids)
+{
     PJ_LP out;
 
-    out.lam = HUGE_VAL; out.phi = HUGE_VAL;
+    out.lam = HUGE_VAL;
+    out.phi = HUGE_VAL;
 
-    ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
-
-    if (ct == nullptr || ct->cvs == nullptr) {
-        if( P->gridlist_count == 1 &&
-            strcmp(P->gridlist[0]->gridname, "null") == 0) {
-          // TODO: remove this particular case that is put there just to be
-          // able to handle longitudes outside of -180,180
-          out = lp;
-        } else {
-          pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        }
+    const auto grid = findGrid(grids, lp);
+    if( !grid ) {
+        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
         return out;
     }
+    if( grid->isNullGrid() )
+    {
+        return lp;
+    }
 
-    inverse = direction == PJ_FWD ? 0 : 1;
-    out = nad_cvt(lp, inverse, ct);
+    int inverse = direction == PJ_FWD ? 0 : 1;
+    out = nad_cvt(lp, inverse, grid);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
-        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
+        pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);
 
     return out;
+}
 
+PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) {
+    return proj_hgrid_apply_internal(P->ctx, lp, direction, P->hgrids);
+}
+
+/************************************************************************/
+/*                         pj_apply_gridshift()                         */
+/*                                                                      */
+/*      This is the externally callable interface - part of the         */
+/*      public API - though it is not used internally any more and I    */
+/*      doubt it is used by any other applications.  But we preserve    */
+/*      it to honour our public api.                                    */
+/************************************************************************/
+
+int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
+                        long point_count, int point_offset,
+                        double *x, double *y, double * /*z */ )
+
+{
+    auto hgrids = getListOfGridSets(ctx, nadgrids);
+    if( hgrids.empty() )
+    {
+        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+        return 1;
+    }
+
+
+    for( long i = 0; i < point_count; i++ )
+    {
+        PJ_LP   input;
+
+        long io = i * point_offset;
+        input.phi = y[io];
+        input.lam = x[io];
+
+        auto output = proj_hgrid_apply_internal(ctx, input, inverse ? PJ_INV : PJ_FWD, hgrids);
+
+        if ( output.lam != HUGE_VAL )
+        {
+            y[io] = output.phi;
+            x[io] = output.lam;
+        }
+        else
+        {
+            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,
+                    y[io] * RAD_TO_DEG );
+            }
+        }
+    }
+
+    return 0;
 }