Skip to content

Commit

Permalink
Merge 73a65df into 366ab09
Browse files Browse the repository at this point in the history
  • Loading branch information
@rouault
rouault committed Dec 24, 2019
2 parents 366ab09 + 73a65df commit 539db9e
Show file tree
Hide file tree
Showing 3 changed files with 73 additions and 38 deletions.
102 changes: 64 additions & 38 deletions src/4D_api.cpp
View file
Expand Up @@ -193,62 +193,88 @@ 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 ");
msg += alt.name;
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
{
Expand All @@ -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 ) {
Expand All @@ -276,7 +303,6 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
}
}
}
i++;
}

proj_errno_set (P, EINVAL);
Expand Down
3 changes: 3 additions & 0 deletions test/cli/ntv2_out.dist
View file
Expand Up @@ -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
6 changes: 6 additions & 0 deletions test/cli/testntv2
View file
Expand Up @@ -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!
Expand Down

0 comments on commit 539db9e

Please sign in to comment.