Add support for chaining vertical CRS transformations with intermediate vertical CRSs sharing a datum with the

source or target

Author: phaarnes

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -560,6 +560,7 @@ struct CoordinateOperationFactory::Private {
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
         bool inCreateOperationsGeogToVertWithAlternativeGeog = false;
         bool inCreateOperationsGeogToVertWithIntermediateVert = false;
+        bool inCreateOperationsVertToVertWithIntermediateVert = false;
         bool skipHorizontalTransformation = false;
         int nRecLevelCreateOperations = 0;
         std::map<std::pair<io::AuthorityFactory::ObjectType, std::string>,
@@ -649,6 +650,15 @@ struct CoordinateOperationFactory::Private {
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
         Context &context);
 
+    static std::vector<CoordinateOperationNNPtr>
+    createOperationsVertToVertWithIntermediateVert(
+        const crs::CRSNNPtr &sourceCRS,
+        const util::optional<common::DataEpoch> &sourceEpoch,
+        const crs::CRSNNPtr &targetCRS,
+        const util::optional<common::DataEpoch> &targetEpoch,
+        const crs::VerticalCRS *vertSrc, const crs::VerticalCRS *vertDst,
+        Context &context);
+
     static void createOperationsFromDatabaseWithVertCRS(
         const crs::CRSNNPtr &sourceCRS,
         const util::optional<common::DataEpoch> &sourceEpoch,
@@ -4530,6 +4540,143 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
 
 // ---------------------------------------------------------------------------
 
+// When transforming between two vertical CRS with different datums, check
+// if there are registered operations whose target (or source) is a different
+// vertical CRS sharing the same datum as our target (or source).  If so,
+// chain: source → intermediate + intermediate → target (the latter being a
+// simple unit/axis conversion handled by createOperationsVertToVert).
+//
+// Typical example: Baltic 1977 height (EPSG:5705) → Caspian depth (EPSG:5706).
+// EPSG has an operation 5705 → 5611 (Caspian height).  5611 and 5706 share the
+// same datum but differ only in axis direction.  We compose:
+//   5705 →(EPSG:5438)→ 5611 →(height-to-depth)→ 5706
+std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
+    createOperationsVertToVertWithIntermediateVert(
+        const crs::CRSNNPtr &sourceCRS,
+        const util::optional<common::DataEpoch> &sourceEpoch,
+        const crs::CRSNNPtr &targetCRS,
+        const util::optional<common::DataEpoch> &targetEpoch,
+        const crs::VerticalCRS *vertSrc, const crs::VerticalCRS *vertDst,
+        Private::Context &context) {
+
+    ENTER_FUNCTION();
+
+    std::vector<CoordinateOperationNNPtr> res;
+
+    struct AntiRecursionGuard {
+        Context &context;
+
+        explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
+            assert(!context.inCreateOperationsVertToVertWithIntermediateVert);
+            context.inCreateOperationsVertToVertWithIntermediateVert = true;
+        }
+
+        ~AntiRecursionGuard() {
+            context.inCreateOperationsVertToVertWithIntermediateVert = false;
+        }
+    };
+    AntiRecursionGuard guard(context);
+
+    const auto &authFactory = context.context->getAuthorityFactory();
+    if (!authFactory)
+        return res;
+    const auto dbContext = authFactory->databaseContext().as_nullable();
+
+    // Strategy 1: pivot through candidates sharing the TARGET's datum.
+    // Find vertical CRSs on the same datum as the target, then look for
+    // operations from the source to each candidate (using the full
+    // createOperations machinery so multi-step chains are also found).
+    auto candidatesDst = findCandidateVertCRSForDatum(
+        authFactory, vertDst->datumNonNull(dbContext).get());
+    for (const auto &candidateVert : candidatesDst) {
+        if (candidateVert->_isEquivalentTo(
+                targetCRS.get(),
+                util::IComparable::Criterion::EQUIVALENT)) {
+            continue; // Skip the target itself — already tried by the caller.
+        }
+        auto opsFirst = createOperations(sourceCRS, sourceEpoch, candidateVert,
+                                         sourceEpoch, context);
+        if (!opsFirst.empty()) {
+            const auto opsSecond = createOperations(
+                candidateVert, sourceEpoch, targetCRS, targetEpoch, context);
+            if (!opsSecond.empty()) {
+                // The transformation from candidateVert to targetCRS should
+                // be just a unit/axis change typically, so take only the
+                // first one, which is likely/hopefully the only one.
+                for (const auto &opFirst : opsFirst) {
+                    if (hasIdentifiers(opFirst)) {
+                        if (candidateVert->_isEquivalentTo(
+                                targetCRS.get(),
+                                util::IComparable::Criterion::EQUIVALENT)) {
+                            res.emplace_back(opFirst);
+                        } else {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opFirst, opsSecond.front()},
+                                            context
+                                                .disallowEmptyIntersection()));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+                if (!res.empty())
+                    return res;
+            }
+        }
+    }
+
+    // Strategy 2: pivot through candidates sharing the SOURCE's datum.
+    // Find vertical CRSs on the same datum as the source, then look for
+    // operations from each candidate to the target.
+    auto candidatesSrc = findCandidateVertCRSForDatum(
+        authFactory, vertSrc->datumNonNull(dbContext).get());
+    for (const auto &candidateVert : candidatesSrc) {
+        if (candidateVert->_isEquivalentTo(
+                sourceCRS.get(),
+                util::IComparable::Criterion::EQUIVALENT)) {
+            continue;
+        }
+        auto opsSecond = createOperations(candidateVert, sourceEpoch, targetCRS,
+                                          targetEpoch, context);
+        if (!opsSecond.empty()) {
+            const auto opsFirst = createOperations(
+                sourceCRS, sourceEpoch, candidateVert, sourceEpoch, context);
+            if (!opsFirst.empty()) {
+                for (const auto &opSecond : opsSecond) {
+                    if (hasIdentifiers(opSecond)) {
+                        if (candidateVert->_isEquivalentTo(
+                                sourceCRS.get(),
+                                util::IComparable::Criterion::EQUIVALENT)) {
+                            res.emplace_back(opSecond);
+                        } else {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opsFirst.front(), opSecond},
+                                            context
+                                                .disallowEmptyIntersection()));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+                if (!res.empty())
+                    return res;
+            }
+        }
+    }
+
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
     createOperationsGeogToVertWithAlternativeGeog(
         const crs::CRSNNPtr &sourceCRS, // geographic CRS
@@ -4663,6 +4810,19 @@ void CoordinateOperationFactory::Private::
         res = applyInverse(res);
     }
 
+    // Typically to transform between two vertical CRS with different datums
+    // (e.g. "Baltic 1977 height" to "Caspian depth") by pivoting through
+    // an intermediate vertical CRS sharing the target's (or source's) datum.
+    // This handles the case where a registered CT targets a height CRS but
+    // the user requests the corresponding depth CRS (or vice versa).
+    if (res.empty() &&
+        !context.inCreateOperationsVertToVertWithIntermediateVert && vertSrc &&
+        vertDst) {
+        res = createOperationsVertToVertWithIntermediateVert(
+            sourceCRS, sourceEpoch, targetCRS, targetEpoch, vertSrc, vertDst,
+            context);
+    }
+
     // There's no direct transformation from NAVD88 height to WGS84,
     // so try to research all transformations from NAVD88 to another
     // intermediate GeographicCRS.

test/unit/test_operationfactory.cpp

@@ -7275,6 +7275,34 @@ TEST(operation, vertCRS_to_vertCRS_New_Zealand_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, vertCRS_to_vertCRS_height_depth_pivot_context) {
+    // Test that PROJ can chain a registered vertical CT with a
+    // height↔depth axis conversion when the target CRS differs from
+    // the CT's registered target only by axis direction.
+    //
+    // EPSG:5705 (Baltic 1977 height) → EPSG:5706 (Caspian depth)
+    // should compose: EPSG:5438 (5705→5611, geogoffset +dh=28)
+    //               + height-to-depth (axisswap order=1,2,-3)
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // Baltic 1977 height
+        authFactory->createCoordinateReferenceSystem("5705"),
+        // Caspian depth (same datum as Caspian height 5611, but axis: down)
+        authFactory->createCoordinateReferenceSystem("5706"), ctxt);
+    ASSERT_GE(list.size(), 1U);
+    // Must NOT be a ballpark transformation
+    EXPECT_FALSE(list[0]->hasBallparkTransformation());
+    // The pipeline should contain geogoffset +dh=28 and axisswap
+    auto projStr =
+        list[0]->exportToPROJString(PROJStringFormatter::create().get());
+    EXPECT_TRUE(projStr.find("geogoffset") != std::string::npos) << projStr;
+    EXPECT_TRUE(projStr.find("axisswap") != std::string::npos) << projStr;
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, projCRS_3D_to_geogCRS_3D) {
 
     auto compoundcrs_ft_obj = PROJStringParser().createFromPROJString(