[Backport 9.8] Chain vertical CRS transformations through intermediate same-datum vertical CRS

docs/source/operations/operations_computation.rst

@@ -566,7 +566,15 @@ Overall logic is similar to the above case. There might be direct operations in
 the PROJ database, involving grid transformations or simple offsets. The fallback
 case is to synthesize a ballpark transformation.
 
-This is implemented by the ``createOperationsVertToVert`` method
+When no direct operation is found, a search is made for intermediate vertical CRS
+entries sharing the same datum as the source or target. If a registered operation
+exists between the source (or target) and such an intermediate CRS, it is composed
+with the axis/unit conversion to the actual target (or from the actual source).
+This is implemented by the ``createOperationsVertToVertWithIntermediateVert``
+method and avoids ballpark fallback for CRS pairs that differ only in axis
+direction (height vs depth), units (e.g. metres vs feet), or both.
+
+The direct case is handled by the ``createOperationsVertToVert`` method.
 
 .. code-block:: shell
 

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>,
@@ -653,6 +654,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,
@@ -4592,15 +4602,17 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
     auto candidatesVert = findCandidateVertCRSForDatum(
         authFactory, vertDst->datumNonNull(dbContext).get());
     for (const auto &candidateVert : candidatesVert) {
+        // Collect all registered source-to-candidate operations, which
+        // may differ in accuracy or area of use, for later ranking.
         auto resTmp = createOperations(sourceCRS, sourceEpoch, candidateVert,
                                        sourceEpoch, context);
         if (!resTmp.empty()) {
+            // The candidate-to-target conversion is expected to be a
+            // trivial unit or axis change, so we only use the first
+            // result (opsSecond.front()).
             const auto opsSecond = createOperations(
                 candidateVert, sourceEpoch, targetCRS, targetEpoch, context);
             if (!opsSecond.empty()) {
-                // The transformation from candidateVert to targetCRS should
-                // be just a unit change typically, so take only the first one,
-                // which is likely/hopefully the only one.
                 for (const auto &opFirst : resTmp) {
                     if (hasIdentifiers(opFirst)) {
                         if (candidateVert->_isEquivalentTo(
@@ -4626,6 +4638,136 @@ 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 to intermediate + intermediate to target (the latter being a
+// simple unit/axis conversion handled by createOperationsVertToVert).
+//
+// Typical example: Baltic 1977 height (EPSG:5705) to Caspian depth (EPSG:5706).
+// EPSG has an operation 5705 to 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();
+
+    // Try both pivot strategies and collect all candidates.  The caller's
+    // filterAndSort will rank them (by accuracy, area of use, etc.)
+
+    // Strategy 1: pivot through candidates sharing the TARGET's datum.
+    // source -> candidate (registered CT) -> target (trivial axis/unit change)
+    {
+        auto candidatesDst = findCandidateVertCRSForDatum(
+            authFactory, vertDst->datumNonNull(dbContext).get());
+        for (const auto &candidateVert : candidatesDst) {
+            if (candidateVert->_isEquivalentTo(
+                    targetCRS.get(),
+                    util::IComparable::Criterion::EQUIVALENT)) {
+                continue;
+            }
+            // Collect all registered source-to-candidate operations, which
+            // may differ in accuracy or area of use, for later ranking.
+            auto opsFirst = createOperations(
+                sourceCRS, sourceEpoch, candidateVert, sourceEpoch, context);
+            if (!opsFirst.empty()) {
+                // The candidate-to-target conversion is expected to be a
+                // trivial unit or axis change, so we only use the first
+                // result (opsSecond.front()).
+                const auto opsSecond =
+                    createOperations(candidateVert, sourceEpoch, targetCRS,
+                                     targetEpoch, context);
+                if (!opsSecond.empty()) {
+                    for (const auto &opFirst : opsFirst) {
+                        if (hasIdentifiers(opFirst)) {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opFirst, opsSecond.front()},
+                                            context
+                                                .disallowEmptyIntersection()));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // 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.
+    if (res.empty()) {
+        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)) {
+                            try {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opsFirst.front(), opSecond},
+                                            context
+                                                .disallowEmptyIntersection()));
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
 std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
     createOperationsGeogToVertWithAlternativeGeog(
         const crs::CRSNNPtr &sourceCRS, // geographic CRS
@@ -4759,6 +4901,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,82 @@ TEST(operation, vertCRS_to_vertCRS_New_Zealand_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, vertCRS_to_vertCRS_pivot_context) {
+    // Test that PROJ can chain a registered vertical CT with a
+    // height-to-depth axis conversion when the target CRS differs from
+    // the CT's registered target only by axis direction.
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+
+    auto checkPipeline = [&](const std::string &srcCode,
+                             const std::string &tgtCode,
+                             const std::string &expectedProj) {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactory->createCoordinateReferenceSystem(srcCode),
+            authFactory->createCoordinateReferenceSystem(tgtCode), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        EXPECT_FALSE(list[0]->hasBallparkTransformation());
+        EXPECT_EQ(
+            list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+            expectedProj);
+    };
+
+    // Using Strategy 1 of createOperationsVertToVertWithIntermediateVert()
+
+    // Caspian: EPSG:5705 (Baltic 1977 height) -> EPSG:5706 (Caspian depth)
+    // via EPSG:5438 (dh=28) + height-to-depth axisswap
+    checkPipeline("5705", "5706",
+                  "+proj=pipeline +step +proj=geogoffset +dh=28 "
+                  "+step +proj=axisswap +order=1,2,-3");
+
+    // Caspian: EPSG:5706 (Caspian depth) -> EPSG:5705 (Baltic 1977 height)
+    // axisswap + inverse of EPSG:5438 (dh=-28)
+    // Note: yes that pipeline is identical to the above one, since it is its
+    // own inverse.
+    checkPipeline("5706", "5705",
+                  "+proj=pipeline +step +proj=geogoffset +dh=28 "
+                  "+step +proj=axisswap +order=1,2,-3");
+
+    // KOC ft: EPSG:5790 -> EPSG:5614 (KOC WD depth ft)
+    // via EPSG:7987 (dh=-4.74) + axisswap + unit conversion m->ft
+    checkPipeline("5790", "5614",
+                  "+proj=pipeline "
+                  "+step +proj=geogoffset +dh=-4.74 "
+                  "+step +proj=axisswap +order=1,2,-3 "
+                  "+step +proj=unitconvert +z_in=m +z_out=ft");
+
+    // KOC ft: EPSG:5614 (KOC WD depth ft) -> EPSG:5790 (KOC CD height)
+    // unit conversion ft->m + axisswap + inverse of EPSG:7987 (dh=4.74)
+    checkPipeline("5614", "5790",
+                  "+proj=pipeline "
+                  "+step +proj=unitconvert +z_in=ft +z_out=m "
+                  "+step +proj=axisswap +order=1,2,-3 "
+                  "+step +proj=geogoffset +dh=4.74");
+
+    // EPSG:5705 (Baltic 1977 height) to EPSG:5336 (Black Sea depth)
+    // Strategy 1 composes: EPSG:5447 (5705 to 5735, geogoffset +dh=0.4)
+    //                    + height-to-depth (axisswap order=1,2,-3)
+    checkPipeline("5705", "5336",
+                  "+proj=pipeline "
+                  "+step +proj=geogoffset +dh=0.4 "
+                  "+step +proj=axisswap +order=1,2,-3");
+
+    // Using Strategy 2 of createOperationsVertToVertWithIntermediateVert()
+
+    // EPSG:5336 (Black Sea depth) to EPSG:5705 (Baltic 1977 height)
+    // Strategy 2 composes: depth-to-height (5336 to 5735, axisswap)
+    //                    + inverse of EPSG:5447 (5735 to 5705, dh=-0.4)
+    checkPipeline("5336", "5705",
+                  "+proj=pipeline "
+                  "+step +proj=axisswap +order=1,2,-3 "
+                  "+step +proj=geogoffset +dh=-0.4");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, projCRS_3D_to_geogCRS_3D) {
 
     auto compoundcrs_ft_obj = PROJStringParser().createFromPROJString(