From f36d2374268608d276da45c96a7b4792ab8ffdb5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 14 Dec 2019 00:16:14 +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/apply_gridshift.cpp | 51 ++++++++++++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 17 deletions(-) (limited to 'src/apply_gridshift.cpp') 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 #include +#include + #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::max(); + dif.phi = std::numeric_limits::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); -- cgit v1.2.3