From 6a4f9cc673e247e1ca76e94d7b1aaaf3e0648c13 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 14 Dec 2019 15:40:17 +0100 Subject: 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 ``` --- src/nad_cvt.cpp | 43 ++++++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 11 deletions(-) (limited to 'src/nad_cvt.cpp') 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 +#include + #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::max(); + dif.phi = std::numeric_limits::max(); + continue; + } + dif.lam = t.lam - del.lam - tb.lam; dif.phi = t.phi + del.phi - tb.phi; t.lam -= dif.lam; -- cgit v1.2.3