diff options
| -rw-r--r-- | src/4D_api.cpp | 102 | ||||
| -rw-r--r-- | test/cli/ntv2_out.dist | 3 | ||||
| -rwxr-xr-x | test/cli/testntv2 | 6 |
3 files changed, 73 insertions, 38 deletions
diff --git a/src/4D_api.cpp b/src/4D_api.cpp index f37594b5..127ca71e 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -193,47 +193,64 @@ similarly, but prefers the 2D resp. 3D interfaces if available. direction = opposite_direction(direction); if( !P->alternativeCoordinateOperations.empty() ) { - // Do a first pass and select the operations that match the area of use - // and has the best accuracy. - int iBest = -1; - double bestAccuracy = std::numeric_limits<double>::max(); - int i = 0; - for( const auto &alt: P->alternativeCoordinateOperations ) { - bool spatialCriterionOK = false; - if( direction == PJ_FWD ) { - if( coord.xyzt.x >= alt.minxSrc && - coord.xyzt.y >= alt.minySrc && - coord.xyzt.x <= alt.maxxSrc && - coord.xyzt.y <= alt.maxySrc) { - spatialCriterionOK = true; + constexpr int N_MAX_RETRY = 2; + int iExcluded[N_MAX_RETRY] = {-1, -1}; + + const int nOperations = static_cast<int>( + P->alternativeCoordinateOperations.size()); + + // We may need several attempts. For example the point at + // lon=-111.5 lat=45.26 falls into the bounding box of the Canadian + // ntv2_0.gsb grid, except that it is not in any of the subgrids, being + // in the US. We thus need another retry that will select the conus + // grid. + for( int iRetry = 0; iRetry <= N_MAX_RETRY; iRetry++ ) + { + // Do a first pass and select the operations that match the area of use + // and has the best accuracy. + int iBest = -1; + double bestAccuracy = std::numeric_limits<double>::max(); + for( int i = 0; i < nOperations; i++ ) { + if( i == iExcluded[0] || i == iExcluded[1] ) { + continue; } - } else { - if( coord.xyzt.x >= alt.minxDst && - coord.xyzt.y >= alt.minyDst && - coord.xyzt.x <= alt.maxxDst && - coord.xyzt.y <= alt.maxyDst ) { - spatialCriterionOK = true; + const auto &alt = P->alternativeCoordinateOperations[i]; + bool spatialCriterionOK = false; + if( direction == PJ_FWD ) { + if( coord.xyzt.x >= alt.minxSrc && + coord.xyzt.y >= alt.minySrc && + coord.xyzt.x <= alt.maxxSrc && + coord.xyzt.y <= alt.maxySrc) { + spatialCriterionOK = true; + } + } else { + if( coord.xyzt.x >= alt.minxDst && + coord.xyzt.y >= alt.minyDst && + coord.xyzt.x <= alt.maxxDst && + coord.xyzt.y <= alt.maxyDst ) { + spatialCriterionOK = true; + } } - } - if( spatialCriterionOK ) { - // The offshore test is for the "Test bug 245 (use +datum=carthage)" - // of testvarious. The long=10 lat=34 point belongs both to the - // onshore and offshore Tunisia area of uses, but is slightly - // onshore. So in a general way, prefer a onshore area to a - // offshore one. - if( iBest < 0 || - (alt.accuracy >= 0 && alt.accuracy < bestAccuracy && - !alt.isOffshore) ) { - iBest = i; - bestAccuracy = alt.accuracy; + if( spatialCriterionOK ) { + // The offshore test is for the "Test bug 245 (use +datum=carthage)" + // of testvarious. The long=10 lat=34 point belongs both to the + // onshore and offshore Tunisia area of uses, but is slightly + // onshore. So in a general way, prefer a onshore area to a + // offshore one. + if( iBest < 0 || + (alt.accuracy >= 0 && alt.accuracy < bestAccuracy && + !alt.isOffshore) ) { + iBest = i; + bestAccuracy = alt.accuracy; + } } } - i ++; - } + if( iBest < 0 ) { + break; + } - if( iBest >= 0 ) { const auto& alt = P->alternativeCoordinateOperations[iBest]; if( P->iCurCoordOp != iBest ) { std::string msg("Using coordinate operation "); @@ -241,14 +258,23 @@ similarly, but prefers the 2D resp. 3D interfaces if available. pj_log(P->ctx, PJ_LOG_TRACE, msg.c_str()); P->iCurCoordOp = iBest; } - return direction == PJ_FWD ? + PJ_COORD res = direction == PJ_FWD ? pj_fwd4d( coord, alt.pj ) : pj_inv4d( coord, alt.pj ); + if( res.xyzt.x != HUGE_VAL ) { + return res; + } + pj_log(P->ctx, PJ_LOG_TRACE, + "Did not result in valid result. " + "Attempting a retry with another operation."); + if( iRetry == N_MAX_RETRY ) { + break; + } + iExcluded[iRetry] = iBest; } // In case we did not find an operation whose area of use is compatible // with the input coordinate, then goes through again the list, and // use the first operation that does not require grids. - i = 0; NS_PROJ::io::DatabaseContextPtr dbContext; try { @@ -257,7 +283,8 @@ similarly, but prefers the 2D resp. 3D interfaces if available. } } catch( const std::exception& ) {} - for( const auto &alt: P->alternativeCoordinateOperations ) { + for( int i = 0; i < nOperations; i++ ) { + const auto &alt = P->alternativeCoordinateOperations[i]; auto coordOperation = dynamic_cast< NS_PROJ::operation::CoordinateOperation*>(alt.pj->iso_obj.get()); if( coordOperation ) { @@ -276,7 +303,6 @@ similarly, but prefers the 2D resp. 3D interfaces if available. } } } - i++; } proj_errno_set (P, EINVAL); diff --git a/test/cli/ntv2_out.dist b/test/cli/ntv2_out.dist index 650a69d8..940997fb 100644 --- a/test/cli/ntv2_out.dist +++ b/test/cli/ntv2_out.dist @@ -12,3 +12,6 @@ Try with NTv2 and NTv1 together ... falls back to NTv1 ############################################################## Switching between NTv2 subgrids -112.5839956 49.4914451 0 -112.58307487 49.49145197 0.00000000 +############################################################## +Attempt first with ntv2_0.gsb and then conus +-111.5 45.26 -111.50079772 45.25992835 0.00000000 diff --git a/test/cli/testntv2 b/test/cli/testntv2 index 2a31304e..ab72d199 100755 --- a/test/cli/testntv2 +++ b/test/cli/testntv2 @@ -60,6 +60,12 @@ $EXE +proj=latlong +datum=NAD83 +to +proj=latlong +ellps=clrk66 +nadgrids=ntv2_0 -112.5839956 49.4914451 0 EOF +echo "##############################################################" >> ${OUT} +echo Attempt first with ntv2_0.gsb and then conus >> ${OUT} +$EXE +proj=longlat +datum=NAD27 +to +proj=longlat +datum=WGS84 -E -d 8 >>${OUT} <<EOF +-111.5 45.26 +EOF + # ############################################################################## # Done! |