createOperations(): fix incorrect height transformation between 3D promoted RGF93 and CH1903+ (fixes #2541)

4 files changed, 122 insertions, 16 deletions

diff --git a/src/iso19111/operation/coordinateoperation_private.hpp b/src/iso19111/operation/coordinateoperation_private.hpp
index 42bedd4b..2dc70414 100644
--- a/src/iso19111/operation/coordinateoperation_private.hpp
+++ b/src/iso19111/operation/coordinateoperation_private.hpp
@@ -50,7 +50,6 @@ struct CoordinateOperation::Private {
     util::optional<common::DataEpoch> sourceCoordinateEpoch_{};
     util::optional<common::DataEpoch> targetCoordinateEpoch_{};
     bool hasBallparkTransformation_ = false;
-    bool use3DHelmert_ = false;
 
     // do not set this for a ProjectedCRS.definingConversion
     struct CRSStrongRef {
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 44d9b570..4832c7b8 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -2614,9 +2614,13 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
                                      bool isNullFirst) {
         const auto opsSecond =
             createOperations(candidateSrcGeod, candidateDstGeod, context);
-        const auto opsThird =
-            createOperations(candidateDstGeod, targetCRS, context);
+        const auto opsThird = createOperations(
+            sourceAndTargetAre3D
+                ? candidateDstGeod->promoteTo3D(std::string(), dbContext)
+                : candidateDstGeod,
+            targetCRS, context);
         assert(!opsThird.empty());
+        const CoordinateOperationNNPtr &opThird(opsThird[0]);
 
         for (auto &opSecond : opsSecond) {
             // Check that it is not a transformation synthetized by
@@ -2632,8 +2636,7 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
             }
 
             std::vector<CoordinateOperationNNPtr> subOps;
-            const bool isNullThird =
-                isNullTransformation(opsThird[0]->nameStr());
+            const bool isNullThird = isNullTransformation(opThird->nameStr());
             CoordinateOperationNNPtr opSecondCloned(
                 (isNullFirst || isNullThird || sourceAndTargetAre3D)
                     ? opSecond->shallowClone()
@@ -2665,12 +2668,31 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
                 }
             }
             if (sourceAndTargetAre3D) {
-                opSecondCloned->getPrivate()->use3DHelmert_ = true;
-                auto invCO = dynamic_cast<InverseCoordinateOperation *>(
-                    opSecondCloned.get());
-                if (invCO) {
-                    auto invCOForward = invCO->forwardOperation().get();
-                    invCOForward->getPrivate()->use3DHelmert_ = true;
+
+                // Force Helmert operations to use the 3D domain, even if the
+                // ones we found in EPSG are advertized for the 2D domain.
+                auto concat =
+                    dynamic_cast<ConcatenatedOperation *>(opSecondCloned.get());
+                if (concat) {
+                    std::vector<CoordinateOperationNNPtr> newSteps;
+                    for (const auto &step : concat->operations()) {
+                        auto newStep = step->shallowClone();
+                        setCRSs(newStep.get(),
+                                newStep->sourceCRS()->promoteTo3D(std::string(),
+                                                                  dbContext),
+                                newStep->targetCRS()->promoteTo3D(std::string(),
+                                                                  dbContext));
+                        newSteps.emplace_back(newStep);
+                    }
+                    opSecondCloned =
+                        ConcatenatedOperation::createComputeMetadata(
+                            newSteps, disallowEmptyIntersection);
+                } else {
+                    setCRSs(opSecondCloned.get(),
+                            opSecondCloned->sourceCRS()->promoteTo3D(
+                                std::string(), dbContext),
+                            opSecondCloned->targetCRS()->promoteTo3D(
+                                std::string(), dbContext));
                 }
             }
             if (isNullFirst) {
@@ -2685,7 +2707,7 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
                 subOps.emplace_back(opSecondCloned);
             } else {
                 subOps.emplace_back(opSecondCloned);
-                subOps.emplace_back(opsThird[0]);
+                subOps.emplace_back(opThird);
             }
 #ifdef TRACE_CREATE_OPERATIONS
             std::string debugStr;
@@ -2726,6 +2748,10 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
     for (const auto &candidateSrcGeod : candidatesSrcGeod) {
         if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
+            auto sourceSrcGeodModified(
+                sourceAndTargetAre3D
+                    ? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
+                    : candidateSrcGeod);
             for (const auto &candidateDstGeod : candidatesDstGeod) {
                 if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
 #ifdef TRACE_CREATE_OPERATIONS
@@ -2734,8 +2760,8 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
                                 objectAsStr(candidateDstGeod.get()) + "->" +
                                 objectAsStr(targetCRS.get()) + ")");
 #endif
-                    const auto opsFirst =
-                        createOperations(sourceCRS, candidateSrcGeod, context);
+                    const auto opsFirst = createOperations(
+                        sourceCRS, sourceSrcGeodModified, context);
                     assert(!opsFirst.empty());
                     const bool isNullFirst =
                         isNullTransformation(opsFirst[0]->nameStr());
@@ -2758,8 +2784,12 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 #ifdef TRACE_CREATE_OPERATIONS
         ENTER_BLOCK("");
 #endif
+        auto sourceSrcGeodModified(
+            sourceAndTargetAre3D
+                ? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
+                : candidateSrcGeod);
         const auto opsFirst =
-            createOperations(sourceCRS, candidateSrcGeod, context);
+            createOperations(sourceCRS, sourceSrcGeodModified, context);
         assert(!opsFirst.empty());
         const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());
 
diff --git a/src/iso19111/operation/transformation.cpp b/src/iso19111/operation/transformation.cpp
index 984f2756..392d6026 100644
--- a/src/iso19111/operation/transformation.cpp
+++ b/src/iso19111/operation/transformation.cpp
@@ -2595,7 +2595,6 @@ void Transformation::_exportToPROJString(
         auto targetCRSGeog =
             dynamic_cast<const crs::GeographicCRS *>(targetCRS().get());
         const bool addPushPopV3 =
-            !CoordinateOperation::getPrivate()->use3DHelmert_ &&
             ((sourceCRSGeog &&
               sourceCRSGeog->coordinateSystem()->axisList().size() == 2) ||
              (targetCRSGeog &&
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index a91d1906..8cc4fa40 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -1453,6 +1453,84 @@ TEST(operation, geocentricCRS_to_geogCRS_different_datum_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_3D_to_geogCRS_3D_different_datum_context) {
+    // Test for https://github.com/OSGeo/PROJ/issues/2541
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // RGF93 (3D)
+        authFactory->createCoordinateReferenceSystem("4965"),
+        // CH1903+ promoted to 3D
+        authFactory->createCoordinateReferenceSystem("4150")->promoteTo3D(
+            std::string(), dbContext),
+        ctxt);
+    ASSERT_GE(list.size(), 1U);
+    EXPECT_EQ(list[0]->nameStr(),
+              "RGF93 to ETRS89 (1) + Inverse of CH1903+ to ETRS89 (1)");
+    // Check that there is no +push +v_3
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+              "+step +proj=cart +ellps=GRS80 "
+              "+step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 "
+              "+step +inv +proj=cart +ellps=bessel "
+              "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+              "+step +proj=axisswap +order=2,1");
+    EXPECT_EQ(list[0]->inverse()->exportToPROJString(
+                  PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+              "+step +proj=cart +ellps=bessel "
+              "+step +proj=helmert +x=674.374 +y=15.056 +z=405.346 "
+              "+step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, geocentric_to_geogCRS_3D_different_datum_context) {
+    // Test variant of https://github.com/OSGeo/PROJ/issues/2541
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // RGF93 (geocentric)
+        authFactory->createCoordinateReferenceSystem("4964"),
+        // CH1903+ promoted to 3D
+        authFactory->createCoordinateReferenceSystem("4150")->promoteTo3D(
+            std::string(), dbContext),
+        ctxt);
+    ASSERT_GE(list.size(), 1U);
+    EXPECT_EQ(list[0]->nameStr(),
+              "Conversion from RGF93 (geocentric) to RGF93 (geog3D) + "
+              "RGF93 to ETRS89 (1) + "
+              "Inverse of CH1903+ to ETRS89 (1)");
+    // Check that there is no +push +v_3
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 "
+              "+step +inv +proj=cart +ellps=bessel "
+              "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+              "+step +proj=axisswap +order=2,1");
+    EXPECT_EQ(list[0]->inverse()->exportToPROJString(
+                  PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
+              "+step +proj=cart +ellps=bessel "
+              "+step +proj=helmert +x=674.374 +y=15.056 +z=405.346");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, esri_projectedCRS_to_geogCRS_with_ITRF_intermediate_context) {
     auto dbContext = DatabaseContext::create();
     auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");