From 9f908ae47cfa70d3cdb2709a8ab5d8eeb10034fc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 9 Dec 2019 23:26:39 +0100 Subject: Horizontal shift grids: hide the 'positive longitude shift value = westward correction' in the CTable2/NTv1/NTv2 backends --- src/apply_gridshift.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'src/apply_gridshift.cpp') diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp index a0ffa394..c786a50a 100644 --- a/src/apply_gridshift.cpp +++ b/src/apply_gridshift.cpp @@ -212,12 +212,12 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) { return t; if (!inverse) { - in.lam -= t.lam; + in.lam += t.lam; in.phi += t.phi; return in; } - t.lam = tb.lam + t.lam; + t.lam = tb.lam - t.lam; t.phi = tb.phi - t.phi; do { @@ -235,8 +235,8 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) { if (del.lam == HUGE_VAL) break; - dif.lam = t.lam - del.lam - tb.lam; - dif.phi = t.phi + del.phi - tb.phi; + dif.lam = t.lam + del.lam - tb.lam; + dif.phi = t.phi + del.phi - tb.phi; t.lam -= dif.lam; t.phi -= dif.phi; -- cgit v1.2.3 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 From 830f94a8117eff270acd2ef928850a9de5e164a9 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 29 Dec 2019 23:57:23 +0100 Subject: proj_hgrid_value(): do not apply compensation for west-oriented longitude_offset. Fixes regression due to grid refactoring work on proj=deformation use case --- src/apply_gridshift.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) (limited to 'src/apply_gridshift.cpp') diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp index b994474b..4ef86fc0 100644 --- a/src/apply_gridshift.cpp +++ b/src/apply_gridshift.cpp @@ -119,7 +119,7 @@ ListOfHGrids proj_hgrid_init(PJ* P, const char *gridkey) { typedef struct { pj_int32 lam, phi; } ILP; -static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid) { +static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid, bool compensateNTConvention) { PJ_LP val, frct; ILP indx; int in; @@ -164,10 +164,10 @@ static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid) { float f10Lon = 0, f10Lat = 0; float f01Lon = 0, f01Lat = 0; float f11Lon = 0, f11Lat = 0; - if( !grid->valueAt(indx.lam, indx.phi, f00Lon, f00Lat) || - !grid->valueAt(indx.lam + 1, indx.phi, f10Lon, f10Lat) || - !grid->valueAt(indx.lam, indx.phi + 1, f01Lon, f01Lat) || - !grid->valueAt(indx.lam + 1, indx.phi + 1, f11Lon, f11Lat) ) + if( !grid->valueAt(indx.lam, indx.phi, compensateNTConvention, f00Lon, f00Lat) || + !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Lon, f10Lat) || + !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Lon, f01Lat) || + !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention, f11Lon, f11Lat) ) { return val; } @@ -210,7 +210,7 @@ PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, const HorizontalShiftGrid* gri tb.lam = adjlon (tb.lam - M_PI) + M_PI; - t = nad_intr (tb, grid); + t = nad_intr (tb, grid, true); if (t.lam == HUGE_VAL) return t; @@ -224,7 +224,7 @@ PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, const HorizontalShiftGrid* gri t.phi = tb.phi - t.phi; do { - del = nad_intr(t, grid); + del = nad_intr(t, grid, true); /* We can possibly go outside of the initial guessed grid, so try */ /* to fetch a new grid into which iterate... */ @@ -297,7 +297,7 @@ PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids& grids, PJ_LP lp) { lp.lam = adjlon(lp.lam - M_PI) + M_PI; - out = nad_intr(lp, grid); + out = nad_intr(lp, grid, false); if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); -- cgit v1.2.3