Skip to content

Commit

Permalink
+proj=gridshift: do not apply iteration in the reverse path with +int…
Browse files Browse the repository at this point in the history 
…erpolation=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.
  • Loading branch information
@rouault
rouault committed Jan 3, 2024
1 parent 5ebe855 commit 6e1294e
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 49 deletions.
118 changes: 70 additions & 48 deletions src/transformations/gridshift.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -549,57 +554,72 @@ 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 easily
// get convergence failures. 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<double>::max();
diff.phi = std::numeric_limits<double>::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<double>::max();
diff.phi = std::numeric_limits<double>::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;
Expand Down Expand Up @@ -667,7 +687,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;
}

Expand Down
15 changes: 14 additions & 1 deletion test/gie/gridshift.gie
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
-------------
Expand Down

0 comments on commit 6e1294e

Please sign in to comment.