createOperations(): when tranforming from a compoundCRS whose vertical component is a BoundCRS, do not apply the horizontal transformation twice

2 files changed, 76 insertions, 19 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d5b766ad..f0ebaebd 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10172,6 +10172,7 @@ struct CoordinateOperationFactory::Private {
         const CoordinateOperationContextNNPtr &context;
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
         bool inCreateOperationsGeogToVertWithIntermediate = false;
+        bool skipHorizontalTransformation = false;
 
         Context(const crs::CRSNNPtr &sourceCRSIn,
                 const crs::CRSNNPtr &targetCRSIn,
@@ -12599,27 +12600,32 @@ CoordinateOperationFactory::Private::createOperations(
             hubSrcGeog->coordinateSystem()->axisList().size() == 3 &&
             geogDst->coordinateSystem()->axisList().size() == 3) {
             auto opsFirst = createOperations(sourceCRS, hubSrc, context);
-            auto opsSecond = createOperations(hubSrc, targetCRS, context);
-            if (!opsFirst.empty() && !opsSecond.empty()) {
-                for (const auto &opFirst : opsFirst) {
-                    for (const auto &opLast : opsSecond) {
-                        // Exclude artificial transformations from the hub
-                        // to the target CRS
-                        if (!opLast->hasBallparkTransformation()) {
-                            try {
-                                res.emplace_back(
-                                    ConcatenatedOperation::
-                                        createComputeMetadata(
-                                            {opFirst, opLast},
-                                            !allowEmptyIntersection));
-                            } catch (
-                                const InvalidOperationEmptyIntersection &) {
+            if (context.skipHorizontalTransformation) {
+                if (!opsFirst.empty())
+                    return opsFirst;
+            } else {
+                auto opsSecond = createOperations(hubSrc, targetCRS, context);
+                if (!opsFirst.empty() && !opsSecond.empty()) {
+                    for (const auto &opFirst : opsFirst) {
+                        for (const auto &opLast : opsSecond) {
+                            // Exclude artificial transformations from the hub
+                            // to the target CRS
+                            if (!opLast->hasBallparkTransformation()) {
+                                try {
+                                    res.emplace_back(
+                                        ConcatenatedOperation::
+                                            createComputeMetadata(
+                                                {opFirst, opLast},
+                                                !allowEmptyIntersection));
+                                } catch (
+                                    const InvalidOperationEmptyIntersection &) {
+                                }
                             }
                         }
                     }
-                }
-                if (!res.empty()) {
-                    return res;
+                    if (!res.empty()) {
+                        return res;
+                    }
                 }
             }
         }
@@ -12841,6 +12847,22 @@ CoordinateOperationFactory::Private::createOperations(
             std::vector<CoordinateOperationNNPtr> verticalTransforms;
             if (componentsSrc.size() >= 2 &&
                 componentsSrc[1]->extractVerticalCRS()) {
+
+                struct SetSkipHorizontalTransform {
+                    Context &context;
+
+                    explicit SetSkipHorizontalTransform(Context &contextIn)
+                        : context(contextIn) {
+                        assert(!context.skipHorizontalTransformation);
+                        context.skipHorizontalTransformation = true;
+                    }
+
+                    ~SetSkipHorizontalTransform() {
+                        context.skipHorizontalTransformation = false;
+                    }
+                };
+                SetSkipHorizontalTransform setSkipHorizontalTransform(context);
+
                 verticalTransforms =
                     createOperations(componentsSrc[1], targetCRS, context);
                 bool foundRegisteredTransformWithAllGridsAvailable = false;
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 88890b10..d52e9bea 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -6014,8 +6014,9 @@ TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) {
         src, NN_NO_CHECK(dst), ctxt);
     ASSERT_EQ(list.size(), 1U);
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+proj=pipeline "
               "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
               "+step +inv +proj=vgridshift +grids=naptrans2008.gtx "
               "+multiplier=1 "
               "+step +inv +proj=hgridshift +grids=rdtrans2008.gsb "
@@ -6026,6 +6027,40 @@ TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, WGS84_to_compoundCRS_with_bound_vertCRS) {
+    auto authFactoryEPSG =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    // WGS 84 3D
+    auto src = authFactoryEPSG->createCoordinateReferenceSystem("4979");
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=NAD83 +geoidgrids=@foo.gtx +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), std::string());
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        src, NN_NO_CHECK(dst), ctxt);
+    ASSERT_GE(list.size(), 53U);
+    EXPECT_EQ(list[0]->nameStr(),
+              "Inverse of unknown to WGS84 ellipsoidal height + "
+              "Inverse of NAD83 to WGS 84 (1) + "
+              "Inverse of axis order change (2D)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +inv +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
 static VerticalCRSNNPtr createVerticalCRS() {
     PropertyMap propertiesVDatum;
     propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG")