Skip to content

Commit

Permalink
proj_crs_to_crs(): tune ordering of operations based on areas, and al…
Browse files Browse the repository at this point in the history 
…so for operations with unknown accuracy (fixes OSGeo#3854)
  • Loading branch information
@rouault
rouault committed Sep 3, 2023
1 parent ed367b8 commit 89943df
Show file tree
Hide file tree
Showing 5 changed files with 47 additions and 26 deletions.
44 changes: 23 additions & 21 deletions src/4D_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,22 +322,14 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
// 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 ||
// If two operations have the same accuracy, use
// the one that is contained within a larger one
(alt.accuracy == bestAccuracy &&
alt.minxSrc >= opList[iBest].minxSrc &&
alt.minySrc >= opList[iBest].minySrc &&
alt.maxxSrc <= opList[iBest].maxxSrc &&
alt.maxySrc <= opList[iBest].maxySrc &&
// check that this is not equality
!(alt.minxSrc == opList[iBest].minxSrc &&
alt.minySrc == opList[iBest].minySrc &&
alt.maxxSrc == opList[iBest].maxxSrc &&
alt.maxySrc == opList[iBest].maxySrc) &&
!opList[iBest].isPriorityOp)) &&
!alt.isOffshore)) {
if (iBest < 0 ||
(((alt.accuracy >= 0 && alt.accuracy < bestAccuracy) ||
// If two operations have the same accuracy, use
// the one that has the smallest area
(alt.accuracy == bestAccuracy &&
alt.pseudoArea < opList[iBest].pseudoArea &&
!opList[iBest].isPriorityOp)) &&
!alt.isOffshore)) {

if (skipNonInstantiable && !alt.isInstantiable()) {
continue;
Expand Down Expand Up @@ -1713,6 +1705,16 @@ static PJ *add_coord_op_to_list(
double maxxDst;
double maxyDst;

double w = west_lon / 180 * M_PI;
double s = south_lat / 180 * M_PI;
double e = east_lon / 180 * M_PI;
double n = north_lat / 180 * M_PI;
if (w > e) {
e += 2 * M_PI;
}
// Integrate cos(lat) between south_lat and north_lat
const double pseudoArea = (e - w) * (std::sin(n) - std::sin(s));

if (pjSrcGeocentricToLonLat) {
minxSrc = west_lon;
minySrc = south_lat;
Expand Down Expand Up @@ -1740,8 +1742,8 @@ static PJ *add_coord_op_to_list(
const double accuracy = proj_coordoperation_get_accuracy(op->ctx, op);
altCoordOps.emplace_back(
idxInOriginalList, minxSrc, minySrc, maxxSrc, maxySrc, minxDst,
minyDst, maxxDst, maxyDst, op, name, accuracy, isOffshore,
pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
minyDst, maxxDst, maxyDst, op, name, accuracy, pseudoArea,
isOffshore, pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
op = nullptr;
}
return op;
Expand Down Expand Up @@ -2987,13 +2989,13 @@ PJCoordOperation::PJCoordOperation(
int idxInOriginalListIn, double minxSrcIn, double minySrcIn,
double maxxSrcIn, double maxySrcIn, double minxDstIn, double minyDstIn,
double maxxDstIn, double maxyDstIn, PJ *pjIn, const std::string &nameIn,
double accuracyIn, bool isOffshoreIn, const PJ *pjSrcGeocentricToLonLatIn,
const PJ *pjDstGeocentricToLonLatIn)
double accuracyIn, double pseudoAreaIn, bool isOffshoreIn,
const PJ *pjSrcGeocentricToLonLatIn, const PJ *pjDstGeocentricToLonLatIn)
: idxInOriginalList(idxInOriginalListIn), minxSrc(minxSrcIn),
minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn),
minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn),
maxyDst(maxyDstIn), pj(pjIn), name(nameIn), accuracy(accuracyIn),
isOffshore(isOffshoreIn),
pseudoArea(pseudoAreaIn), isOffshore(isOffshoreIn),
isPriorityOp(isSpecialCaseForNAD83_to_NAD83HARN(name) ||
isSpecialCaseForGDA94_to_WGS84(name) ||
isSpecialCaseForWGS84_to_GDA2020(name)),
Expand Down
3 changes: 2 additions & 1 deletion src/iso19111/c_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9110,7 +9110,8 @@ PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) {
alt.idxInOriginalList, minxSrc, minySrc, maxxSrc,
maxySrc, minxDst, minyDst, maxxDst, maxyDst,
pjNormalized, co->nameStr(), alt.accuracy,
alt.isOffshore, alt.pjSrcGeocentricToLonLat,
alt.pseudoArea, alt.isOffshore,
alt.pjSrcGeocentricToLonLat,
alt.pjDstGeocentricToLonLat);
}
}
Expand Down
10 changes: 6 additions & 4 deletions src/proj_internal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,6 +342,7 @@ struct PJCoordOperation {
PJ *pj = nullptr;
std::string name{};
double accuracy = -1.0;
double pseudoArea = 0.0;
bool isOffshore = false;
bool isPriorityOp = false;
bool srcIsLonLatDegree = false;
Expand All @@ -363,7 +364,7 @@ struct PJCoordOperation {
double minySrcIn, double maxxSrcIn, double maxySrcIn,
double minxDstIn, double minyDstIn, double maxxDstIn,
double maxyDstIn, PJ *pjIn, const std::string &nameIn,
double accuracyIn, bool isOffshoreIn,
double accuracyIn, double pseudoAreaIn, bool isOffshoreIn,
const PJ *pjSrcGeocentricToLonLatIn,
const PJ *pjDstGeocentricToLonLatIn);

Expand All @@ -376,7 +377,8 @@ struct PJCoordOperation {
minyDst(other.minyDst), maxxDst(other.maxxDst),
maxyDst(other.maxyDst), pj(proj_clone(ctx, other.pj)),
name(std::move(other.name)), accuracy(other.accuracy),
isOffshore(other.isOffshore), isPriorityOp(other.isPriorityOp),
pseudoArea(other.pseudoArea), isOffshore(other.isOffshore),
isPriorityOp(other.isPriorityOp),
srcIsLonLatDegree(other.srcIsLonLatDegree),
srcIsLatLonDegree(other.srcIsLatLonDegree),
dstIsLonLatDegree(other.dstIsLonLatDegree),
Expand All @@ -396,8 +398,8 @@ struct PJCoordOperation {
maxySrc(other.maxySrc), minxDst(other.minxDst),
minyDst(other.minyDst), maxxDst(other.maxxDst),
maxyDst(other.maxyDst), name(std::move(other.name)),
accuracy(other.accuracy), isOffshore(other.isOffshore),
isPriorityOp(other.isPriorityOp),
accuracy(other.accuracy), pseudoArea(other.pseudoArea),
isOffshore(other.isOffshore), isPriorityOp(other.isPriorityOp),
srcIsLonLatDegree(other.srcIsLonLatDegree),
srcIsLatLonDegree(other.srcIsLatLonDegree),
dstIsLonLatDegree(other.dstIsLonLatDegree),
Expand Down
12 changes: 12 additions & 0 deletions test/cli/testvarious
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1320,6 +1320,18 @@ $EXE -d 3 EPSG:25831 EPSG:23031 -E >>${OUT} 2>&1 <<EOF
299905.060 4499796.515
EOF

echo "##############################################################" >> ${OUT}
echo "Test Similarity Transformation through CompoundCRS" >> ${OUT}
# Cf https://github.com/OSGeo/PROJ/issues/3854#issuecomment-1689964773
#
$EXE -d 3 EPSG:3912 EPSG:3794 -E >>${OUT} 2>&1 <<EOF
477134.28 95134.21
EOF

$EXE -d 3 EPSG:3912+EPSG:5779 EPSG:3794+EPSG:8690 -E >>${OUT} 2>&1 <<EOF
477134.28 95134.21 5
EOF


# Done!
# do 'diff' with distribution results
Expand Down
4 changes: 4 additions & 0 deletions test/cli/tv_out.dist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -633,3 +633,7 @@ Test Similarity Transformation (example from EPSG Guidance Note 7.2)
##############################################################
Test inverse of Similarity Transformation (example from EPSG Guidance Note 7.2)
299905.060 4499796.515 300000.000 4500000.000 0.000
##############################################################
Test Similarity Transformation through CompoundCRS
477134.28 95134.21 476763.303 95620.222 0.000
477134.28 95134.21 5 476763.303 95620.222 5.000

0 comments on commit 89943df

Please sign in to comment.