From 5ebe855f43bb74474d128ef5132cdbb09d7ae220 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 3 Jan 2024 14:14:45 +0100 Subject: [PATCH 1/2] proj.h: add PROJ_ERR_COORD_TRANSFM_NO_CONVERGENCE error code --- src/proj.h | 2 ++ src/strerrno.cpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/proj.h b/src/proj.h index cef28e3ccc..f0543ddfca 100644 --- a/src/proj.h +++ b/src/proj.h @@ -706,6 +706,8 @@ PJ_COORD PROJ_DLL proj_geod(const PJ *P, PJ_COORD a, PJ_COORD b); #define PROJ_ERR_COORD_TRANSFM_GRID_AT_NODATA \ (PROJ_ERR_COORD_TRANSFM + \ 5) /* point to transform falls in a grid cell that evaluates to nodata */ +#define PROJ_ERR_COORD_TRANSFM_NO_CONVERGENCE \ + (PROJ_ERR_COORD_TRANSFM + 6) /* iterative convergence method fail */ /** Other type of errors */ #define PROJ_ERR_OTHER 4096 diff --git a/src/strerrno.cpp b/src/strerrno.cpp index 2e35935d43..0784f5b6b5 100644 --- a/src/strerrno.cpp +++ b/src/strerrno.cpp @@ -34,6 +34,8 @@ static const struct { {PROJ_ERR_COORD_TRANSFM_GRID_AT_NODATA, _("Coordinate to transform falls into a grid cell that evaluates to " "nodata")}, + {PROJ_ERR_COORD_TRANSFM_NO_CONVERGENCE, + _("Iterative method fails to converge on coordinate to transform")}, {PROJ_ERR_OTHER_API_MISUSE, _("API misuse")}, {PROJ_ERR_OTHER_NO_INVERSE_OP, _("No inverse operation")}, {PROJ_ERR_OTHER_NETWORK_ERROR, From c7d41bb876887e03aaace61a13f4d5d467b02160 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 3 Jan 2024 14:15:41 +0100 Subject: [PATCH 2/2] +proj=gridshift: do not apply iteration in the reverse path with +interpolation=biquadratic As NOAA NCAT transformation tool does. Otherwise we may get convergence errors on points that are close to the boundary of cells or half-cells. --- ..._nad83_harn_conus_extract_sanfrancisco.tif | Bin 0 -> 1285 bytes src/transformations/gridshift.cpp | 119 +++++++++++------- test/gie/gridshift.gie | 15 ++- 3 files changed, 85 insertions(+), 49 deletions(-) create mode 100644 data/tests/us_noaa_nadcon5_nad83_1986_nad83_harn_conus_extract_sanfrancisco.tif diff --git a/data/tests/us_noaa_nadcon5_nad83_1986_nad83_harn_conus_extract_sanfrancisco.tif b/data/tests/us_noaa_nadcon5_nad83_1986_nad83_harn_conus_extract_sanfrancisco.tif new file mode 100644 index 0000000000000000000000000000000000000000..c2acae8a2e0de62ce3649dd198ccfda27b3a14f4 GIT binary patch literal 1285 zcmZ`&TWlLe6rFWK2*FfTwD18%qeY+ugm>+v!QM&SxN+LJO>HYS67a!l>>WF+*t=$T z?LJf>l!t$KilP=&_()5oN~LWOs=&ibR6;>WRP=)a;!z=#50EGisS*zf?ra7Qtf?_g4P#Z+{1Z)Nb6d5c$Z!nzrxVHO z-~je)yiA5SM%rZmWMVo!Fc8Ayw$)&cxAL_a&#SqbqLi6!Fn4#zvMCLf?S@jT7Zi&c zMcdk_h_9+KHSDkoqb^f6*UgIIir3LWJsgV;hnC?Cb*vCX)Ln!|wk8u(JD5ie>d|-~ z#v>_@Rk1~@Y*eNWqq<$v>n@YycjprLPLY?#oT_Qrb$3J&+W!ZcLQu0S)H7{MuQG4O zHspB0e4tJZ2gqXQFpoNA=IKSp)mhCg>qS_o93O(M>8#J+o1Kv3lbPw%eVKGVF{Nj= zY{^aJJByPOx$)^#HlNC*W$e;wtpcl4W$f4jz>XFjvnKMD<05abZZKGP$z@(=znblu z9@K)~pp`oN*D@@}VdpOs)G>*0=Ok>a+$oC@<-flRp_g713&{`1qfT-AH3DPpenXM}$G{m1}dX;gkR!&(L>vhlXR5x?7MH0i`^-Ha{EVQ+tMuWzvff2R^7#)+9X`tTkfWi){I6G^BY%H>fFC#J$@tF&UJEXe4VMr%E?q}%{`w4`eIiKy&`)#r z=T7qMfv5O`kDMdh#t%gs=g*MXGr88LIbHkVyCwepiI?AcUVDQtK75Fr*!6nr?4I+a zpI2HR-uf}QOI^pO=Eulq7y915asQXv+_%5-FTPqN?3WK(%Grw~b;nd|edz-EMC;-I E02UZ$C;$Ke literal 0 HcmV?d00001 diff --git a/src/transformations/gridshift.cpp b/src/transformations/gridshift.cpp index bda2e3d415..44504ae30d 100644 --- a/src/transformations/gridshift.cpp +++ b/src/transformations/gridshift.cpp @@ -80,7 +80,8 @@ struct gridshiftData { const PJ_LPZ &input, GenericShiftGridSet *&gridSetOut) const; PJ_LPZ grid_interpolate(PJ_CONTEXT *ctx, const std::string &type, PJ_LP lp, - const GenericShiftGrid *grid); + const GenericShiftGrid *grid, + bool &biquadraticInterpolationOut); PJ_LPZ grid_apply_internal(PJ_CONTEXT *ctx, const std::string &type, bool isVerticalOnly, const PJ_LPZ in, PJ_DIRECTION direction, @@ -178,7 +179,8 @@ typedef struct { #define REL_TOLERANCE_HGRIDSHIFT 1e-5 PJ_LPZ gridshiftData::grid_interpolate(PJ_CONTEXT *ctx, const std::string &type, - PJ_LP lp, const GenericShiftGrid *grid) { + PJ_LP lp, const GenericShiftGrid *grid, + bool &biquadraticInterpolationOut) { PJ_LPZ val; val.lam = val.phi = HUGE_VAL; @@ -274,6 +276,7 @@ PJ_LPZ gridshiftData::grid_interpolate(PJ_CONTEXT *ctx, const std::string &type, const int idxSampleLong = gridInfo.idxSampleLong; const int idxSampleZ = gridInfo.idxSampleZ; const bool bilinearInterpolation = gridInfo.bilinearInterpolation; + biquadraticInterpolationOut = !bilinearInterpolation; ILP indx; const auto &extent = grid->extentAndRes(); @@ -521,7 +524,9 @@ PJ_LPZ gridshiftData::grid_apply_internal( const NS_PROJ::ExtentAndRes *extent; PJ_LP normalized_in = normalizeLongitude(grid, in, extent); - PJ_LPZ shift = grid_interpolate(ctx, type, normalized_in, grid); + bool biquadraticInterpolationOut = false; + PJ_LPZ shift = grid_interpolate(ctx, type, normalized_in, grid, + biquadraticInterpolationOut); if (grid->hasChanged()) { shouldRetry = gridset->reopen(ctx); PJ_LPZ out; @@ -549,57 +554,73 @@ PJ_LPZ gridshiftData::grid_apply_internal( guess.lam = normalized_in.lam - shift.lam; guess.phi = normalized_in.phi - shift.phi; - int i = MAX_ITERATIONS; - const double toltol = TOL * TOL; - PJ_LP diff; - do { - shift = grid_interpolate(ctx, type, guess, grid); - if (grid->hasChanged()) { - shouldRetry = gridset->reopen(ctx); + // NOAA NCAT transformer tool doesn't do iteration in the reverse path. + // Do the same (only for biquadratic, although NCAT applies this logic to + // bilinear too) + // Cf + // https://github.com/noaa-ngs/ncat-lib/blob/77bcff1ce4a78fe06d0312102ada008aefcc2c62/src/gov/noaa/ngs/grid/Transformer.java#L374 + // When trying to do iterative reverse path with biquadratic, we can + // get convergence failures on points that are close to the boundary of + // cells or half-cells. For example with + // echo -122.4250009683 37.8286740788 0 | bin/cct -I +proj=gridshift + // +grids=tests/us_noaa_nadcon5_nad83_1986_nad83_harn_conus_extract_sanfrancisco.tif + // +interpolation=biquadratic + if (!biquadraticInterpolationOut) { + int i = MAX_ITERATIONS; + const double toltol = TOL * TOL; + PJ_LP diff; + do { + shift = grid_interpolate(ctx, type, guess, grid, + biquadraticInterpolationOut); + if (grid->hasChanged()) { + shouldRetry = gridset->reopen(ctx); + PJ_LPZ out; + out.lam = out.phi = out.z = HUGE_VAL; + return out; + } + + /* We can possibly go outside of the initial guessed grid, so try */ + /* to fetch a new grid into which iterate... */ + if (shift.lam == HUGE_VAL) { + PJ_LPZ lp; + lp.lam = guess.lam; + lp.phi = guess.phi; + auto newGrid = findGrid(type, lp, gridset); + if (newGrid == nullptr || newGrid == grid || + newGrid->isNullGrid()) + break; + pj_log(ctx, PJ_LOG_TRACE, "Switching from grid %s to grid %s", + grid->name().c_str(), newGrid->name().c_str()); + grid = newGrid; + normalized_in = normalizeLongitude(grid, in, extent); + diff.lam = std::numeric_limits::max(); + diff.phi = std::numeric_limits::max(); + continue; + } + + diff.lam = guess.lam + shift.lam - normalized_in.lam; + diff.phi = guess.phi + shift.phi - normalized_in.phi; + guess.lam -= diff.lam; + guess.phi -= diff.phi; + + } while (--i && (diff.lam * diff.lam + diff.phi * diff.phi > + toltol)); /* prob. slightly faster than hypot() */ + + if (i == 0) { + pj_log(ctx, PJ_LOG_TRACE, + "Inverse grid shift iterator failed to converge."); + proj_context_errno_set(ctx, PROJ_ERR_COORD_TRANSFM_NO_CONVERGENCE); PJ_LPZ out; out.lam = out.phi = out.z = HUGE_VAL; return out; } - /* We can possibly go outside of the initial guessed grid, so try */ - /* to fetch a new grid into which iterate... */ if (shift.lam == HUGE_VAL) { - PJ_LPZ lp; - lp.lam = guess.lam; - lp.phi = guess.phi; - auto newGrid = findGrid(type, lp, gridset); - if (newGrid == nullptr || newGrid == grid || newGrid->isNullGrid()) - break; - pj_log(ctx, PJ_LOG_TRACE, "Switching from grid %s to grid %s", - grid->name().c_str(), newGrid->name().c_str()); - grid = newGrid; - normalized_in = normalizeLongitude(grid, in, extent); - diff.lam = std::numeric_limits::max(); - diff.phi = std::numeric_limits::max(); - continue; + pj_log( + ctx, PJ_LOG_TRACE, + "Inverse grid shift iteration failed, presumably at grid edge. " + "Using first approximation."); } - - diff.lam = guess.lam + shift.lam - normalized_in.lam; - diff.phi = guess.phi + shift.phi - normalized_in.phi; - guess.lam -= diff.lam; - guess.phi -= diff.phi; - - } while (--i && (diff.lam * diff.lam + diff.phi * diff.phi > - toltol)); /* prob. slightly faster than hypot() */ - - if (i == 0) { - pj_log(ctx, PJ_LOG_TRACE, - "Inverse grid shift iterator failed to converge."); - proj_context_errno_set(ctx, PROJ_ERR_COORD_TRANSFM); - PJ_LPZ out; - out.lam = out.phi = out.z = HUGE_VAL; - return out; - } - - if (shift.lam == HUGE_VAL) { - pj_log(ctx, PJ_LOG_TRACE, - "Inverse grid shift iteration failed, presumably at grid edge. " - "Using first approximation."); } PJ_LPZ out; @@ -667,7 +688,9 @@ PJ_LPZ gridshiftData::apply(PJ *P, PJ_DIRECTION direction, PJ_LPZ lpz) { } if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { - proj_context_errno_set(P->ctx, PROJ_ERR_COORD_TRANSFM_OUTSIDE_GRID); + if (proj_context_errno(P->ctx) == 0) { + proj_context_errno_set(P->ctx, PROJ_ERR_COORD_TRANSFM_OUTSIDE_GRID); + } return out; } diff --git a/test/gie/gridshift.gie b/test/gie/gridshift.gie index 86de12346e..61d5729b1d 100644 --- a/test/gie/gridshift.gie +++ b/test/gie/gridshift.gie @@ -209,7 +209,7 @@ roundtrip 1 ------------------------------------------------------------------------------- operation +proj=gridshift +grids=tests/us_noaa_nadcon5_nad83_2007_nad83_2011_conus_extract.tif +interpolation=biquadratic ------------------------------------------------------------------------------- -tolerance 0.001 mm +tolerance 0.005 mm accept -95.4916666666 37.0083333333 10.000 expect -95.49166648893 37.00833334837 9.984340 @@ -224,6 +224,19 @@ accept -95.4916666666 37.0083333333 10.000 expect -95.49166648893 37.00833334838 9.984341 roundtrip 1 +---------------------------------------------------------------------------------------------------------------------- +# Test case with inverse biquadratic convergence where we are around a location where the interpolation window changes +---------------------------------------------------------------------------------------------------------------------- + +------------------------------------------------------------------------------- +operation +proj=gridshift +grids=tests/us_noaa_nadcon5_nad83_1986_nad83_harn_conus_extract_sanfrancisco.tif +interpolation=biquadratic +------------------------------------------------------------------------------- +direction inverse +tolerance 0.005 mm + +accept -122.4250009683 37.8286740788 +expect -122.4249999391 37.8286728006 + ------------- # Error cases -------------