Horizontal grid shift: fix issue on iterative inverse computation when switching between (sub)grids (fixes #1663)

Given in.txt with
53.999759140 5.144478208 252.6995

Before the fix,
cct  -t 0 -d 4 +proj=pipeline +step +proj=axisswap +order=2,1,3,4 +step +proj=hgridshift +inv +grids=rdtrans2018.gsb +step +proj=vgridshift +grids=naptrans2018.gtx +step +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel in.txt

returned:
  139079.8814    668306.0302      212.1724        0.0000

It now returns:
  139079.8850    668306.0458      212.1724        0.0000

which meets with the 1mm accuracy the expected result of test point
```
30010049    53.999759140    5.144478208     252.6995        139079.8850     668306.0460     212.1723
```

6 files changed, 51 insertions, 17 deletions

diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
index e6c7b2b1..39e7c3b5 100644
--- a/src/apply_gridshift.cpp
+++ b/src/apply_gridshift.cpp
@@ -115,7 +115,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
 /*    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) {
+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 */
@@ -210,7 +210,7 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_coun
         ct = find_ctable(ctx, input, gridlist_count, gridlist);
         if( ct != nullptr )
         {
-            output = nad_cvt( input, inverse, ct );
+            output = nad_cvt( ctx, input, inverse, ct, gridlist_count, gridlist);
 
             if ( output.lam != HUGE_VAL && debug_count++ < 20 )
                 pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id );
@@ -356,7 +356,7 @@ PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) {
     }
 
     inverse = direction == PJ_FWD ? 0 : 1;
-    out = nad_cvt(lp, inverse, ct);
+    out = nad_cvt(P->ctx, lp, inverse, ct, P->gridlist_count, P->gridlist);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
         pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
diff --git a/src/gridcatalog.cpp b/src/gridcatalog.cpp
index 15d81dd7..9b94fef8 100644
--- a/src/gridcatalog.cpp
+++ b/src/gridcatalog.cpp
@@ -164,7 +164,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse,
             return PJD_ERR_FAILED_TO_LOAD_GRID;
         }
             
-        output_after = nad_cvt( input, inverse, gi->ct );
+        output_after = nad_cvt( defn->ctx, input, inverse, gi->ct, 0, nullptr );
         if( output_after.lam == HUGE_VAL )
         {
             if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
@@ -213,7 +213,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse,
             return PJD_ERR_FAILED_TO_LOAD_GRID;
         }
             
-        output_before = nad_cvt( input, inverse, gi->ct );
+        output_before = nad_cvt( defn->ctx, input, inverse, gi->ct, 0, nullptr );
         if( output_before.lam == HUGE_VAL )
         {
             if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
diff --git a/src/nad_cvt.cpp b/src/nad_cvt.cpp
index 79441d0a..e8b8e9b7 100644
--- a/src/nad_cvt.cpp
+++ b/src/nad_cvt.cpp
@@ -7,10 +7,12 @@
 #include "proj_internal.h"
 #include <math.h>
 
+#include <limits>
+
 #define MAX_ITERATIONS 10
 #define TOL 1e-12
 
-PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) {
+PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables) {
     PJ_LP t, tb,del, dif;
     int i = MAX_ITERATIONS;
     const double toltol = TOL*TOL;
@@ -40,18 +42,37 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) {
     do {
         del = nad_intr(t, ct);
 
-        /* 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)
+        /* In case we are (or have switched) on the null grid, exit now */
+        if( del.lam == 0 && del.phi == 0 )
             break;
 
+        /* We can possibly go outside of the initial guessed grid, so try */
+        /* to fetch a new grid into which iterate... */
+        if (del.lam == HUGE_VAL)
+        {
+            if( grid_count == 0 )
+                break;
+            PJ_LP lp;
+            lp.lam = t.lam + ct->ll.lam;
+            lp.phi = t.phi + ct->ll.phi;
+            auto newCt = find_ctable(ctx, lp, grid_count, tables);
+            if( newCt == nullptr || newCt == ct )
+                break;
+            pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s",
+                   ct->id,
+                   newCt->id);
+            ct = newCt;
+            t.lam = lp.lam - ct->ll.lam;
+            t.phi = lp.phi - ct->ll.phi;
+            tb = in;
+            tb.lam -= ct->ll.lam;
+            tb.phi -= ct->ll.phi;
+            tb.lam = adjlon (tb.lam - M_PI) + M_PI;
+            dif.lam = std::numeric_limits<double>::max();
+            dif.phi = std::numeric_limits<double>::max();
+            continue;
+        }
+
         dif.lam = t.lam - del.lam - tb.lam;
            dif.phi = t.phi + del.phi - tb.phi;
         t.lam -= dif.lam;
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 3ca927a3..3e219682 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -863,8 +863,10 @@ int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *);
 int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *);
 
 /* nadcon related protos */
+struct CTABLE* find_ctable(projCtx_t *ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables);
+
 PJ_LP             nad_intr(PJ_LP, struct CTABLE *);
-PJ_LP             nad_cvt(PJ_LP, int, struct CTABLE *);
+PJ_LP             nad_cvt(projCtx_t *ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables);
 struct CTABLE *nad_init(projCtx_t *ctx, char *);
 struct CTABLE *nad_ctable_init( projCtx_t *ctx, struct projFileAPI_t* fid );
 int            nad_ctable_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid );
diff --git a/test/cli/ntv2_out.dist b/test/cli/ntv2_out.dist
index 9a97f9cf..650a69d8 100644
--- a/test/cli/ntv2_out.dist
+++ b/test/cli/ntv2_out.dist
@@ -9,3 +9,6 @@ Try with NTv2 and NTv1 together ... falls back to NTv1
 99d00'00.000"W 65d00'00.000"N 0.0	99d0'1.5885"W	65d0'1.3482"N 0.000
 111d00'00.000"W 46d00'00.000"N 0.0	111d0'3.1897"W	45d59'59.7489"N 0.000
 111d00'00.000"W 47d30'00.000"N 0.0	111d0'2.7989"W	47d29'59.9896"N 0.000
+##############################################################
+Switching between NTv2 subgrids
+-112.5839956 49.4914451 0	-112.58307487	49.49145197 0.00000000
diff --git a/test/cli/testntv2 b/test/cli/testntv2
index 73371dbe..2a31304e 100755
--- a/test/cli/testntv2
+++ b/test/cli/testntv2
@@ -52,6 +52,14 @@ $EXE +proj=latlong +ellps=clrk66 +nadgrids=ntv2_0.gsb,ntv1_can.dat,conus \
 111d00'00.000"W 46d00'00.000"N 0.0
 111d00'00.000"W 47d30'00.000"N 0.0
 EOF
+
+echo "##############################################################" >> ${OUT}
+echo Switching between NTv2 subgrids >> ${OUT}
+# Initial guess is in ALraymnd, going to parent CAwest afterwards
+$EXE +proj=latlong +datum=NAD83 +to +proj=latlong +ellps=clrk66 +nadgrids=ntv2_0.gsb -E -d 8 >>${OUT} <<EOF
+-112.5839956 49.4914451 0
+EOF
+
 #
 ##############################################################################
 # Done!