Merge pull request #2227 from rouault/fix_2225

createOperations(): when converting CompoundCRS<-->Geographic3DCrs, do not use discard change of ellipsoidal height if a Helmert transformation is involved

1 files changed, 61 insertions, 45 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index ad456860..35be415f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -119,8 +119,6 @@ static const char *BALLPARK_GEOCENTRIC_TRANSLATION =
 static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
 static const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
 static const char *BALLPARK_GEOGRAPHIC_OFFSET = "Ballpark geographic offset";
-static const char *BALLPARK_VERTICAL_TRANSFORMATION_PREFIX =
-    " (ballpark vertical transformation";
 static const char *BALLPARK_VERTICAL_TRANSFORMATION =
     " (ballpark vertical transformation)";
 static const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
@@ -128,6 +126,8 @@ static const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
     "height correction)";
 
 static const std::string AXIS_ORDER_CHANGE_2D_NAME = "axis order change (2D)";
+static const std::string AXIS_ORDER_CHANGE_3D_NAME =
+    "axis order change (geographic3D horizontal)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -4842,8 +4842,7 @@ Conversion::createHeightDepthReversal(const util::PropertyMap &properties) {
  */
 ConversionNNPtr Conversion::createAxisOrderReversal(bool is3D) {
     if (is3D) {
-        return create(createMapNameEPSGCode(
-                          "axis order change (geographic3D horizontal)", 15499),
+        return create(createMapNameEPSGCode(AXIS_ORDER_CHANGE_3D_NAME, 15499),
                       createMethodMapNameEPSGCode(
                           EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_3D),
                       {}, {});
@@ -7829,7 +7828,8 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
                 }
                 if (starts_with(tokens[i], INVERSE_OF)) {
                     name += tokens[i].substr(INVERSE_OF.size());
-                } else if (tokens[i] == AXIS_ORDER_CHANGE_2D_NAME) {
+                } else if (tokens[i] == AXIS_ORDER_CHANGE_2D_NAME ||
+                           tokens[i] == AXIS_ORDER_CHANGE_3D_NAME) {
                     name += tokens[i];
                 } else {
                     name += INVERSE_OF + tokens[i];
@@ -11123,8 +11123,8 @@ struct SortFunction {
 
     // Sorting function
     // Return true if a < b
-    bool operator()(const CoordinateOperationNNPtr &a,
-                    const CoordinateOperationNNPtr &b) const {
+    bool compare(const CoordinateOperationNNPtr &a,
+                 const CoordinateOperationNNPtr &b) const {
         auto iterA = map.find(a.get());
         assert(iterA != map.end());
         auto iterB = map.find(b.get());
@@ -11249,6 +11249,15 @@ struct SortFunction {
         // Arbitrary final criterion
         return a_name < b_name;
     }
+
+    bool operator()(const CoordinateOperationNNPtr &a,
+                    const CoordinateOperationNNPtr &b) const {
+        const bool ret = compare(a, b);
+#if 0
+        std::cerr << a->nameStr() << " < " << b->nameStr() << " : " << ret << std::endl;
+#endif
+        return ret;
+    }
 };
 
 // ---------------------------------------------------------------------------
@@ -11264,6 +11273,17 @@ static size_t getStepCount(const CoordinateOperationNNPtr &op) {
 
 // ---------------------------------------------------------------------------
 
+static bool isNullTransformation(const std::string &name) {
+    if (name.find(" + ") != std::string::npos)
+        return false;
+    return starts_with(name, BALLPARK_GEOCENTRIC_TRANSLATION) ||
+           starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET) ||
+           starts_with(name, NULL_GEOGRAPHIC_OFFSET) ||
+           starts_with(name, NULL_GEOCENTRIC_TRANSLATION);
+}
+
+// ---------------------------------------------------------------------------
+
 struct FilterResults {
 
     FilterResults(const std::vector<CoordinateOperationNNPtr> &sourceListIn,
@@ -11514,19 +11534,6 @@ struct FilterResults {
 
             const auto stepCount = getStepCount(op);
 
-            const bool isApprox =
-                op->nameStr().find(BALLPARK_VERTICAL_TRANSFORMATION_PREFIX) !=
-                std::string::npos;
-            const bool isNullTransformation =
-                op->nameStr().find(BALLPARK_GEOGRAPHIC_OFFSET) !=
-                    std::string::npos ||
-                op->nameStr().find(NULL_GEOGRAPHIC_OFFSET) !=
-                    std::string::npos ||
-                op->nameStr().find(NULL_GEOCENTRIC_TRANSLATION) !=
-                    std::string::npos ||
-                op->nameStr().find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
-                    std::string::npos;
-
             bool isPROJExportable = false;
             auto formatter = io::PROJStringFormatter::create();
             try {
@@ -11537,10 +11544,24 @@ struct FilterResults {
             } catch (const std::exception &) {
             }
 
+#if 0
+            std::cerr << op->nameStr() << " ";
+            std::cerr << area << " ";
+            std::cerr << getAccuracy(op) << " ";
+            std::cerr << isPROJExportable << " ";
+            std::cerr << hasGrids << " ";
+            std::cerr << gridsAvailable << " ";
+            std::cerr << gridsKnown << " ";
+            std::cerr << stepCount << " ";
+            std::cerr << op->hasBallparkTransformation() << " ";
+            std::cerr << isNullTransformation(op->nameStr()) << " ";
+            std::cerr << std::endl;
+#endif
             map[op.get()] = PrecomputedOpCharacteristics(
                 area, getAccuracy(op), isPROJExportable, hasGrids,
-                gridsAvailable, gridsKnown, stepCount, isApprox,
-                isNullTransformation);
+                gridsAvailable, gridsKnown, stepCount,
+                op->hasBallparkTransformation(),
+                isNullTransformation(op->nameStr()));
         }
 
         // Sort !
@@ -11587,12 +11608,8 @@ struct FilterResults {
         // better
         if (hasOpThatContainsAreaOfInterestAndNoGrid && res.size() > 1) {
             const auto &opLast = res.back();
-            const std::string &name = opLast->nameStr();
-            if (name.find(BALLPARK_GEOGRAPHIC_OFFSET) != std::string::npos ||
-                name.find(NULL_GEOGRAPHIC_OFFSET) != std::string::npos ||
-                name.find(NULL_GEOCENTRIC_TRANSLATION) != std::string::npos ||
-                name.find(BALLPARK_GEOCENTRIC_TRANSLATION) !=
-                    std::string::npos) {
+            if (opLast->hasBallparkTransformation() ||
+                isNullTransformation(opLast->nameStr())) {
                 std::vector<CoordinateOperationNNPtr> resTemp;
                 for (size_t i = 0; i < res.size() - 1; i++) {
                     resTemp.emplace_back(res[i]);
@@ -12438,7 +12455,8 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
         horizTransform, verticalTransform, geogDst);
 
     const bool horizTransformIsNoOp =
-        starts_with(horizTransform->nameStr(), NULL_GEOGRAPHIC_OFFSET);
+        starts_with(horizTransform->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+        horizTransform->nameStr().find(" + ") == std::string::npos;
     if (horizTransformIsNoOp) {
         auto properties = util::PropertyMap();
         properties.set(common::IdentifiedObject::NAME_KEY,
@@ -12496,11 +12514,13 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
             interpolationGeogCRS);
 
     std::vector<CoordinateOperationNNPtr> ops;
-    if (!starts_with(opSrcCRSToGeogCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+    if (!(starts_with(opSrcCRSToGeogCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+          opSrcCRSToGeogCRS->nameStr().find(" + ") == std::string::npos)) {
         ops.emplace_back(opSrcCRSToGeogCRS);
     }
     ops.emplace_back(verticalTransform);
-    if (!starts_with(opGeogCRStoDstCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+    if (!(starts_with(opGeogCRStoDstCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+          opGeogCRStoDstCRS->nameStr().find(" + ") == std::string::npos)) {
         ops.emplace_back(opGeogCRStoDstCRS);
     }
 
@@ -12792,16 +12812,6 @@ findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
 
 // ---------------------------------------------------------------------------
 
-static bool isNullTransformation(const std::string &name) {
-
-    return starts_with(name, BALLPARK_GEOCENTRIC_TRANSLATION) ||
-           starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET) ||
-           starts_with(name, NULL_GEOGRAPHIC_OFFSET) ||
-           starts_with(name, NULL_GEOCENTRIC_TRANSLATION);
-}
-
-// ---------------------------------------------------------------------------
-
 void CoordinateOperationFactory::Private::setCRSs(
     CoordinateOperation *co, const crs::CRSNNPtr &sourceCRS,
     const crs::CRSNNPtr &targetCRS) {
@@ -14789,10 +14799,16 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
                         if (!componentsSrc[0]->isEquivalentTo(
                                 target2D.get(),
                                 util::IComparable::Criterion::EQUIVALENT)) {
-                            interpToTargetOps = createOperations(
-                                NN_NO_CHECK(interpolationGeogCRS),
-                                targetCRS->demoteTo2D(std::string(), dbContext),
-                                context);
+                            // We do the transformation from the
+                            // interpolationCRS
+                            // to the target one in 3D (see #2225)
+                            // But we don't do that between sourceCRS and
+                            // interpolationCRS, as this would mess with an
+                            // orthometric elevation.
+                            auto interp3D = interpolationGeogCRS->promoteTo3D(
+                                std::string(), dbContext);
+                            interpToTargetOps =
+                                createOperations(interp3D, targetCRS, context);
                         }
                     };