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/apply_gridshift.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src/apply_gridshift.cpp') 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); -- cgit v1.2.3