diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/apply_gridshift.cpp | 6 | ||||
| -rw-r--r-- | src/gridcatalog.cpp | 4 | ||||
| -rw-r--r-- | src/nad_cvt.cpp | 43 | ||||
| -rw-r--r-- | src/proj_internal.h | 4 |
4 files changed, 40 insertions, 17 deletions
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); diff --git a/src/gridcatalog.cpp b/src/gridcatalog.cpp index 15d81dd7..9b94fef8 100644 --- a/src/gridcatalog.cpp +++ b/src/gridcatalog.cpp @@ -164,7 +164,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse, return PJD_ERR_FAILED_TO_LOAD_GRID; } - output_after = nad_cvt( input, inverse, gi->ct ); + output_after = nad_cvt( defn->ctx, input, inverse, gi->ct, 0, nullptr ); if( output_after.lam == HUGE_VAL ) { if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) @@ -213,7 +213,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse, return PJD_ERR_FAILED_TO_LOAD_GRID; } - output_before = nad_cvt( input, inverse, gi->ct ); + output_before = nad_cvt( defn->ctx, input, inverse, gi->ct, 0, nullptr ); if( output_before.lam == HUGE_VAL ) { if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) diff --git a/src/nad_cvt.cpp b/src/nad_cvt.cpp index 79441d0a..e8b8e9b7 100644 --- a/src/nad_cvt.cpp +++ b/src/nad_cvt.cpp @@ -7,10 +7,12 @@ #include "proj_internal.h" #include <math.h> +#include <limits> + #define MAX_ITERATIONS 10 #define TOL 1e-12 -PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) { +PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables) { PJ_LP t, tb,del, dif; int i = MAX_ITERATIONS; const double toltol = TOL*TOL; @@ -40,18 +42,37 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) { do { del = nad_intr(t, ct); - /* 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. */ - if (del.lam == HUGE_VAL) + /* In case we are (or have switched) on the null grid, exit now */ + if( del.lam == 0 && del.phi == 0 ) break; + /* 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) + { + if( grid_count == 0 ) + break; + PJ_LP lp; + lp.lam = t.lam + ct->ll.lam; + lp.phi = t.phi + ct->ll.phi; + auto newCt = find_ctable(ctx, lp, grid_count, tables); + if( newCt == nullptr || newCt == ct ) + break; + pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s", + ct->id, + newCt->id); + ct = newCt; + t.lam = lp.lam - ct->ll.lam; + t.phi = lp.phi - ct->ll.phi; + tb = in; + tb.lam -= ct->ll.lam; + tb.phi -= ct->ll.phi; + tb.lam = adjlon (tb.lam - M_PI) + M_PI; + dif.lam = std::numeric_limits<double>::max(); + dif.phi = std::numeric_limits<double>::max(); + continue; + } + dif.lam = t.lam - del.lam - tb.lam; dif.phi = t.phi + del.phi - tb.phi; t.lam -= dif.lam; diff --git a/src/proj_internal.h b/src/proj_internal.h index 3ca927a3..3e219682 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -863,8 +863,10 @@ int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *); int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *); /* nadcon related protos */ +struct CTABLE* find_ctable(projCtx_t *ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables); + PJ_LP nad_intr(PJ_LP, struct CTABLE *); -PJ_LP nad_cvt(PJ_LP, int, struct CTABLE *); +PJ_LP nad_cvt(projCtx_t *ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables); struct CTABLE *nad_init(projCtx_t *ctx, char *); struct CTABLE *nad_ctable_init( projCtx_t *ctx, struct projFileAPI_t* fid ); int nad_ctable_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid ); |