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
```

1 files changed, 34 insertions, 17 deletions

diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
index c786a50a..b994474b 100644
--- a/src/apply_gridshift.cpp
+++ b/src/apply_gridshift.cpp
@@ -36,6 +36,8 @@
 #include <stdio.h>
 #include <string.h>
 
+#include <limits>
+
 #include "proj.h"
 #include "proj_internal.h"
 #include "proj/internal/internal.hpp"
@@ -191,7 +193,8 @@ static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid) {
 #define TOL 1e-12
 
 static
-PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) {
+PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, const HorizontalShiftGrid* grid,
+              const ListOfHGrids& grids) {
     PJ_LP t, tb,del, dif;
     int i = MAX_ITERATIONS;
     const double toltol = TOL*TOL;
@@ -201,9 +204,9 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) {
 
     /* normalize input to ll origin */
     tb = in;
-    const auto& extent = grid->extentAndRes();
-    tb.lam -= extent.westLon;
-    tb.phi -= extent.southLat;
+    const auto* extent = &(grid->extentAndRes());
+    tb.lam -= extent->westLon;
+    tb.phi -= extent->southLat;
 
     tb.lam = adjlon (tb.lam - M_PI) + M_PI;
 
@@ -223,17 +226,31 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) {
     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. */
+        /* 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)
-            break;
+        {
+            PJ_LP lp;
+            lp.lam = t.lam + extent->westLon;
+            lp.phi = t.phi + extent->southLat;
+            auto newGrid = findGrid(grids, lp);
+            if( newGrid == nullptr || newGrid == grid || newGrid->isNullGrid() )
+                break;
+            pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s",
+                   grid->name().c_str(),
+                   newGrid->name().c_str());
+            grid = newGrid;
+            extent = &(grid->extentAndRes());
+            t.lam = lp.lam - extent->westLon;
+            t.phi = lp.phi - extent->southLat;
+            tb = in;
+            tb.lam -= extent->westLon;
+            tb.phi -= extent->southLat;
+            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;
@@ -254,8 +271,8 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) {
     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;
+    in.lam = adjlon (t.lam + extent->westLon);
+    in.phi = t.phi + extent->southLat;
     return in;
 }
 
@@ -311,7 +328,7 @@ PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx,
     }
 
     int inverse = direction == PJ_FWD ? 0 : 1;
-    out = nad_cvt(lp, inverse, grid);
+    out = nad_cvt(ctx, lp, inverse, grid, grids);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
         pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);