pj_get_suggested_operation(): tune it to give correct result for RGAF…

…09 to RRAF 1991 / UTM zone 20N + Guadeloupe 1988 height transformation
rouault · Nov 16, 2023 · fbc49aa · fbc49aa
1 parent 3e33c27
commit fbc49aa
Show file tree

Hide file tree

Showing 6 changed files with 36 additions and 15 deletions.
diff --git a/data/tests/fr_ign_RAGTBT2016.tif b/data/tests/fr_ign_RAGTBT2016.tif
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
@@ -328,6 +328,8 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
                   // the one that has the smallest area
                   (alt.accuracy == bestAccuracy &&
                    alt.pseudoArea < opList[iBest].pseudoArea &&
+                   !(alt.isUnknownAreaName &&
+                     !opList[iBest].isUnknownAreaName) &&
                    !opList[iBest].isPriorityOp)) &&
                  !alt.isOffshore)) {
 
@@ -1693,7 +1695,7 @@ static PJ *add_coord_op_to_list(
     int idxInOriginalList, PJ *op, double west_lon, double south_lat,
     double east_lon, double north_lat, PJ *pjGeogToSrc, PJ *pjGeogToDst,
     const PJ *pjSrcGeocentricToLonLat, const PJ *pjDstGeocentricToLonLat,
-    bool isOffshore, std::vector<PJCoordOperation> &altCoordOps) {
+    const char *areaName, std::vector<PJCoordOperation> &altCoordOps) {
     /*****************************************************************************/
 
     double minxSrc;
@@ -1742,8 +1744,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, pseudoArea,
-            isOffshore, pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
+            minyDst, maxxDst, maxyDst, op, name, accuracy, pseudoArea, areaName,
+            pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
         op = nullptr;
     }
     return op;
@@ -2018,23 +2020,22 @@ pj_create_prepared_operations(PJ_CONTEXT *ctx, const PJ *source_crs,
                 north_lat = 90;
             }
 
-            const bool isOffshore = areaName && strstr(areaName, "- offshore");
             if (west_lon <= east_lon) {
                 op = add_coord_op_to_list(
                     i, op, west_lon, south_lat, east_lon, north_lat,
                     pjGeogToSrc, pjGeogToDst, pjSrcGeocentricToLonLat,
-                    pjDstGeocentricToLonLat, isOffshore, preparedOpList);
+                    pjDstGeocentricToLonLat, areaName, preparedOpList);
             } else {
                 auto op_clone = proj_clone(ctx, op);
 
                 op = add_coord_op_to_list(
                     i, op, west_lon, south_lat, 180, north_lat, pjGeogToSrc,
                     pjGeogToDst, pjSrcGeocentricToLonLat,
-                    pjDstGeocentricToLonLat, isOffshore, preparedOpList);
+                    pjDstGeocentricToLonLat, areaName, preparedOpList);
                 op_clone = add_coord_op_to_list(
                     i, op_clone, -180, south_lat, east_lon, north_lat,
                     pjGeogToSrc, pjGeogToDst, pjSrcGeocentricToLonLat,
-                    pjDstGeocentricToLonLat, isOffshore, preparedOpList);
+                    pjDstGeocentricToLonLat, areaName, preparedOpList);
                 proj_destroy(op_clone);
             }
 
@@ -3019,13 +3020,15 @@ 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, double pseudoAreaIn, bool isOffshoreIn,
+    double accuracyIn, double pseudoAreaIn, const char *areaNameIn,
     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),
-      pseudoArea(pseudoAreaIn), isOffshore(isOffshoreIn),
+      pseudoArea(pseudoAreaIn), areaName(areaNameIn ? areaNameIn : ""),
+      isOffshore(areaName.find("- offshore") != std::string::npos),
+      isUnknownAreaName(areaName.empty() || areaName == "unknown"),
       isPriorityOp(isSpecialCaseForNAD83_to_NAD83HARN(name) ||
                    isSpecialCaseForGDA94_to_WGS84(name) ||
                    isSpecialCaseForWGS84_to_GDA2020(name)),

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
@@ -9135,7 +9135,7 @@ 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.pseudoArea, alt.isOffshore,
+                        alt.pseudoArea, alt.areaName.c_str(),
                         alt.pjSrcGeocentricToLonLat,
                         alt.pjDstGeocentricToLonLat);
                 }

diff --git a/src/proj_internal.h b/src/proj_internal.h
@@ -343,7 +343,9 @@ struct PJCoordOperation {
     std::string name{};
     double accuracy = -1.0;
     double pseudoArea = 0.0;
+    std::string areaName{};
     bool isOffshore = false;
+    bool isUnknownAreaName = false;
     bool isPriorityOp = false;
     bool srcIsLonLatDegree = false;
     bool srcIsLatLonDegree = false;
@@ -364,8 +366,8 @@ struct PJCoordOperation {
                      double minySrcIn, double maxxSrcIn, double maxySrcIn,
                      double minxDstIn, double minyDstIn, double maxxDstIn,
                      double maxyDstIn, PJ *pjIn, const std::string &nameIn,
-                     double accuracyIn, double pseudoAreaIn, bool isOffshoreIn,
-                     const PJ *pjSrcGeocentricToLonLatIn,
+                     double accuracyIn, double pseudoAreaIn,
+                     const char *areaName, const PJ *pjSrcGeocentricToLonLatIn,
                      const PJ *pjDstGeocentricToLonLatIn);
 
     PJCoordOperation(const PJCoordOperation &) = delete;
@@ -377,7 +379,9 @@ 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),
-          pseudoArea(other.pseudoArea), isOffshore(other.isOffshore),
+          pseudoArea(other.pseudoArea), areaName(other.areaName),
+          isOffshore(other.isOffshore),
+          isUnknownAreaName(other.isUnknownAreaName),
           isPriorityOp(other.isPriorityOp),
           srcIsLonLatDegree(other.srcIsLonLatDegree),
           srcIsLatLonDegree(other.srcIsLatLonDegree),
@@ -399,7 +403,9 @@ struct PJCoordOperation {
           minyDst(other.minyDst), maxxDst(other.maxxDst),
           maxyDst(other.maxyDst), name(std::move(other.name)),
           accuracy(other.accuracy), pseudoArea(other.pseudoArea),
-          isOffshore(other.isOffshore), isPriorityOp(other.isPriorityOp),
+          areaName(std::move(other.areaName)), isOffshore(other.isOffshore),
+          isUnknownAreaName(other.isUnknownAreaName),
+          isPriorityOp(other.isPriorityOp),
           srcIsLonLatDegree(other.srcIsLonLatDegree),
           srcIsLatLonDegree(other.srcIsLatLonDegree),
           dstIsLonLatDegree(other.dstIsLonLatDegree),
@@ -422,7 +428,7 @@ struct PJCoordOperation {
                maxxDst == other.maxxDst && maxyDst == other.maxyDst &&
                name == other.name &&
                proj_is_equivalent_to(pj, other.pj, PJ_COMP_STRICT) &&
-               accuracy == other.accuracy && isOffshore == other.isOffshore;
+               accuracy == other.accuracy && areaName == other.areaName;
     }
 
     bool operator!=(const PJCoordOperation &other) const {

diff --git a/test/cli/testvarious b/test/cli/testvarious
@@ -1359,6 +1359,14 @@ $EXE -d 3 +proj=topocentric +datum=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 +Z_0=3
 0 0 0
 EOF
 
+echo "##############################################################" >> ${OUT}
+echo "Test cs2cs EPSG:5488 (RGAF09) to EPSG:4559+5757 (RRAF 1991 / UTM zone 20N + Guadeloupe 1988 height)" >> ${OUT}
+echo "Check that we use the horizontal transformation for Guadeloupe (RRAF 1991 to RGAF09 (2)), and not the one for Martinique" >> ${OUT}
+#
+PROJ_DATA=${PROJ_DATA}:${PROJ_DATA}/tests $EXE -d 3 EPSG:5488 EPSG:4559+5757 >> ${OUT} <<EOF
+16.248285304 -61.484212843 53.073
+EOF
+
 
 # Done!
 # do 'diff' with distribution results

diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
@@ -649,3 +649,7 @@ Test --t_epoch
 ##############################################################
 Test +proj=topocentric +datum=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 +Z_0=3692463.58 +no_defs +type=crs +to EPSG:4978
 0 0 0	-3982059.420	3331314.880 3692463.580
+##############################################################
+Test cs2cs EPSG:5488 (RGAF09) to EPSG:4559+5757 (RRAF 1991 / UTM zone 20N + Guadeloupe 1988 height)
+Check that we use the horizontal transformation for Guadeloupe (RRAF 1991 to RGAF09 (2)), and not the one for Martinique
+661991.318	1796999.201 93.846