Vertical transformations: improve situations similar to transforming from 'NAVD88 (ftUS)' to X, where we now consider the available transformations from 'NAVD88' to X that might exist in the database

4 files changed, 465 insertions, 163 deletions

diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index 12b3b111..37b94c1e 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -1089,6 +1089,10 @@ class PROJ_GCC_DLL AuthorityFactory {
                                const std::string &datum_code,
                                const std::string &geodetic_crs_type) const;
 
+    PROJ_INTERNAL std::list<crs::VerticalCRSNNPtr>
+    createVerticalCRSFromDatum(const std::string &datum_auth_name,
+                               const std::string &datum_code) const;
+
     PROJ_INTERNAL std::list<crs::GeodeticCRSNNPtr>
     createGeodeticCRSFromEllipsoid(const std::string &ellipsoid_auth_name,
                                    const std::string &ellipsoid_code,
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index c7581642..27ac712c 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10298,6 +10298,7 @@ struct CoordinateOperationFactory::Private {
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
         bool inCreateOperationsThroughPreferredHub = false;
         bool inCreateOperationsGeogToVertWithAlternativeGeog = false;
+        bool inCreateOperationsGeogToVertWithIntermediateVert = false;
         bool skipHorizontalTransformation = false;
 
         Context(const crs::CRSNNPtr &sourceCRSIn,
@@ -10327,6 +10328,23 @@ struct CoordinateOperationFactory::Private {
         const crs::VerticalCRS *vertDst,
         std::vector<CoordinateOperationNNPtr> &res);
 
+    static std::vector<CoordinateOperationNNPtr>
+    createOperationsGeogToVertWithIntermediateVert(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        const crs::VerticalCRS *vertDst, Context &context);
+
+    static std::vector<CoordinateOperationNNPtr>
+    createOperationsGeogToVertWithAlternativeGeog(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Context &context);
+
+    static void createOperationsFromDatabaseWithVertCRS(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::GeographicCRS *geogSrc,
+        const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
+        const crs::VerticalCRS *vertDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
     static void createOperationsGeodToGeod(
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
         Private::Context &context, const crs::GeodeticCRS *geodSrc,
@@ -10409,11 +10427,6 @@ struct CoordinateOperationFactory::Private {
         const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst,
         Context &context);
 
-    static std::vector<CoordinateOperationNNPtr>
-    createOperationsGeogToVertWithAlternativeGeog(
-        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-        Context &context);
-
     static bool
     hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> &res,
                              const Context &context);
@@ -12402,51 +12415,6 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub(
 
 // ---------------------------------------------------------------------------
 
-std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
-    createOperationsGeogToVertWithAlternativeGeog(
-        const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS
-        const crs::CRSNNPtr &targetCRS,      // vertical CRS
-        Private::Context &context) {
-
-    std::vector<CoordinateOperationNNPtr> res;
-
-    struct AntiRecursionGuard {
-        Context &context;
-
-        explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
-            assert(!context.inCreateOperationsGeogToVertWithAlternativeGeog);
-            context.inCreateOperationsGeogToVertWithAlternativeGeog = true;
-        }
-
-        ~AntiRecursionGuard() {
-            context.inCreateOperationsGeogToVertWithAlternativeGeog = false;
-        }
-    };
-    AntiRecursionGuard guard(context);
-
-    for (int i = 0; i < 2; i++) {
-
-        // Generally EPSG has operations from GeogCrs to VertCRS
-        auto ops =
-            i == 0 ? findOpsInRegistryDirectTo(targetCRS, context.context)
-                   : findOpsInRegistryDirectFrom(targetCRS, context.context);
-
-        for (const auto &op : ops) {
-            const auto tmpCRS = i == 0 ? op->sourceCRS() : op->targetCRS();
-            if (tmpCRS &&
-                dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) {
-                res.emplace_back(i == 0 ? op : op->inverse());
-            }
-        }
-        if (!res.empty())
-            break;
-    }
-
-    return res;
-}
-
-// ---------------------------------------------------------------------------
-
 static CoordinateOperationNNPtr
 createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
                                     const crs::CRSNNPtr &targetCRS) {
@@ -12723,88 +12691,32 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
         doFilterAndCheckPerfectOp = false;
 
         bool sameGeodeticDatum = false;
-        if (geodSrc && geodDst) {
+
+        if (vertSrc || vertDst) {
+            createOperationsFromDatabaseWithVertCRS(sourceCRS, targetCRS,
+                                                    context, geogSrc, geogDst,
+                                                    vertSrc, vertDst, res);
+        } else if (geodSrc && geodDst) {
+
             const auto &srcDatum = geodSrc->datum();
             const auto &dstDatum = geodDst->datum();
             sameGeodeticDatum =
                 srcDatum != nullptr && dstDatum != nullptr &&
                 srcDatum->_isEquivalentTo(
                     dstDatum.get(), util::IComparable::Criterion::EQUIVALENT);
-        }
 
-        // NAD83 only exists in 2D version in EPSG, so if it has been
-        // promotted to 3D, when researching a vertical to geog
-        // transformation, try to down cast to 2D.
-        const auto geog3DToVertTryThroughGeog2D = [&res, &context](
-            const crs::GeographicCRS *geogSrcIn,
-            const crs::VerticalCRS *vertDstIn,
-            const crs::CRSNNPtr &targetCRSIn) {
-            if (res.empty() && geogSrcIn && vertDstIn &&
-                geogSrcIn->coordinateSystem()->axisList().size() == 3 &&
-                geogSrcIn->datum()) {
-                const auto &authFactory =
-                    context.context->getAuthorityFactory();
-                const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
-                    authFactory, geogSrcIn->datum().get()));
-                for (const auto &candidate : candidatesSrcGeod) {
-                    auto geogCandidate =
-                        util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
-                            candidate);
-                    if (geogCandidate &&
-                        geogCandidate->coordinateSystem()->axisList().size() ==
-                            2) {
-                        res = findOpsInRegistryDirect(
-                            NN_NO_CHECK(geogCandidate), targetCRSIn,
-                            context.context);
-                        if (res.empty()) {
-                            res = applyInverse(findOpsInRegistryDirect(
-                                targetCRSIn, NN_NO_CHECK(geogCandidate),
-                                context.context));
-                        }
-                        break;
-                    }
-                }
-                return true;
-            }
-            return false;
-        };
-
-        if (geog3DToVertTryThroughGeog2D(geogSrc, vertDst, targetCRS)) {
-            // do nothing
-        } else if (geog3DToVertTryThroughGeog2D(geogDst, vertSrc, sourceCRS)) {
-            res = applyInverse(res);
-        }
-
-        // There's no direct transformation from NAVD88 height to WGS84,
-        // so try to research all transformations from NAVD88 to another
-        // intermediate GeographicCRS.
-        if (res.empty() &&
-            !context.inCreateOperationsGeogToVertWithAlternativeGeog &&
-            geogSrc && vertDst) {
-            res = createOperationsGeogToVertWithAlternativeGeog(
-                sourceCRS, targetCRS, context);
-        } else if (res.empty() &&
-                   !context.inCreateOperationsGeogToVertWithAlternativeGeog &&
-                   geogDst && vertSrc) {
-            res = applyInverse(createOperationsGeogToVertWithAlternativeGeog(
-                targetCRS, sourceCRS, context));
-        }
-
-        if (res.empty() && !sameGeodeticDatum &&
-            !context.inCreateOperationsWithDatumPivotAntiRecursion && geodSrc &&
-            geodDst) {
-            // If we still didn't find a transformation, and that the source
-            // and target are GeodeticCRS, then go through their underlying
-            // datum to find potential transformations between other
-            // GeodeticRSs
-            // that are made of those datum
-            // The typical example is if transforming between two
-            // GeographicCRS,
-            // but transformations are only available between their
-            // corresponding geocentric CRS.
-            const auto &srcDatum = geodSrc->datum();
-            const auto &dstDatum = geodDst->datum();
-            if (srcDatum != nullptr && dstDatum != nullptr) {
+            if (res.empty() && !sameGeodeticDatum &&
+                !context.inCreateOperationsWithDatumPivotAntiRecursion &&
+                srcDatum != nullptr && dstDatum != nullptr) {
+                // If we still didn't find a transformation, and that the source
+                // and target are GeodeticCRS, then go through their underlying
+                // datum to find potential transformations between other
+                // GeodeticRSs
+                // that are made of those datum
+                // The typical example is if transforming between two
+                // GeographicCRS,
+                // but transformations are only available between their
+                // corresponding geocentric CRS.
                 createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
                                                geodSrc, geodDst, context);
                 doFilterAndCheckPerfectOp = !res.empty();
@@ -12851,6 +12763,222 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
 
 // ---------------------------------------------------------------------------
 
+static std::vector<crs::CRSNNPtr>
+findCandidateVertCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
+                             const datum::VerticalReferenceFrame *datum) {
+    std::vector<crs::CRSNNPtr> candidates;
+    assert(datum);
+    const auto &ids = datum->identifiers();
+    const auto &datumName = datum->nameStr();
+    if (!ids.empty()) {
+        for (const auto &id : ids) {
+            const auto &authName = *(id->codeSpace());
+            const auto &code = id->code();
+            if (!authName.empty()) {
+                auto l_candidates =
+                    authFactory->createVerticalCRSFromDatum(authName, code);
+                for (const auto &candidate : l_candidates) {
+                    candidates.emplace_back(candidate);
+                }
+            }
+        }
+    } else if (datumName != "unknown" && datumName != "unnamed") {
+        auto matches = authFactory->createObjectsFromName(
+            datumName,
+            {io::AuthorityFactory::ObjectType::VERTICAL_REFERENCE_FRAME}, false,
+            2);
+        if (matches.size() == 1) {
+            const auto &match = matches.front();
+            if (datum->_isEquivalentTo(
+                    match.get(), util::IComparable::Criterion::EQUIVALENT) &&
+                !match->identifiers().empty()) {
+                return findCandidateVertCRSForDatum(
+                    authFactory,
+                    dynamic_cast<const datum::VerticalReferenceFrame *>(
+                        match.get()));
+            }
+        }
+    }
+    return candidates;
+}
+
+// ---------------------------------------------------------------------------
+
+std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
+    createOperationsGeogToVertWithIntermediateVert(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        const crs::VerticalCRS *vertDst, Private::Context &context) {
+
+    std::vector<CoordinateOperationNNPtr> res;
+
+    struct AntiRecursionGuard {
+        Context &context;
+
+        explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
+            assert(!context.inCreateOperationsGeogToVertWithIntermediateVert);
+            context.inCreateOperationsGeogToVertWithIntermediateVert = true;
+        }
+
+        ~AntiRecursionGuard() {
+            context.inCreateOperationsGeogToVertWithIntermediateVert = false;
+        }
+    };
+    AntiRecursionGuard guard(context);
+    const auto &authFactory = context.context->getAuthorityFactory();
+    auto candidatesVert =
+        findCandidateVertCRSForDatum(authFactory, vertDst->datum().get());
+    for (const auto &candidateVert : candidatesVert) {
+        auto resTmp = createOperations(sourceCRS, candidateVert, context);
+        if (!resTmp.empty()) {
+            const auto opsSecond =
+                createOperations(candidateVert, targetCRS, 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)) {
+                        res.emplace_back(
+                            ConcatenatedOperation::createComputeMetadata(
+                                {opFirst, opsSecond.front()},
+                                !allowEmptyIntersection));
+                    }
+                }
+                if (!res.empty())
+                    break;
+            }
+        }
+    }
+
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
+std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
+    createOperationsGeogToVertWithAlternativeGeog(
+        const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS
+        const crs::CRSNNPtr &targetCRS,      // vertical CRS
+        Private::Context &context) {
+
+    std::vector<CoordinateOperationNNPtr> res;
+
+    struct AntiRecursionGuard {
+        Context &context;
+
+        explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
+            assert(!context.inCreateOperationsGeogToVertWithAlternativeGeog);
+            context.inCreateOperationsGeogToVertWithAlternativeGeog = true;
+        }
+
+        ~AntiRecursionGuard() {
+            context.inCreateOperationsGeogToVertWithAlternativeGeog = false;
+        }
+    };
+    AntiRecursionGuard guard(context);
+
+    for (int i = 0; i < 2; i++) {
+
+        // Generally EPSG has operations from GeogCrs to VertCRS
+        auto ops =
+            i == 0 ? findOpsInRegistryDirectTo(targetCRS, context.context)
+                   : findOpsInRegistryDirectFrom(targetCRS, context.context);
+
+        for (const auto &op : ops) {
+            const auto tmpCRS = i == 0 ? op->sourceCRS() : op->targetCRS();
+            if (tmpCRS &&
+                dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) {
+                res.emplace_back(i == 0 ? op : op->inverse());
+            }
+        }
+        if (!res.empty())
+            break;
+    }
+
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::
+    createOperationsFromDatabaseWithVertCRS(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::GeographicCRS *geogSrc,
+        const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
+        const crs::VerticalCRS *vertDst,
+        std::vector<CoordinateOperationNNPtr> &res) {
+
+    // Typically to transform from "NAVD88 height (ftUS)" to a geog CRS
+    // by using transformations of "NAVD88 height" (metre) to that geog CRS
+    if (res.empty() &&
+        !context.inCreateOperationsGeogToVertWithIntermediateVert && geogSrc &&
+        vertDst) {
+        res = createOperationsGeogToVertWithIntermediateVert(
+            sourceCRS, targetCRS, vertDst, context);
+    } else if (res.empty() &&
+               !context.inCreateOperationsGeogToVertWithIntermediateVert &&
+               geogDst && vertSrc) {
+        res = applyInverse(createOperationsGeogToVertWithIntermediateVert(
+            targetCRS, sourceCRS, vertSrc, context));
+    }
+
+    // NAD83 only exists in 2D version in EPSG, so if it has been
+    // promotted to 3D, when researching a vertical to geog
+    // transformation, try to down cast to 2D.
+    const auto geog3DToVertTryThroughGeog2D = [&res, &context](
+        const crs::GeographicCRS *geogSrcIn, const crs::VerticalCRS *vertDstIn,
+        const crs::CRSNNPtr &targetCRSIn) {
+        if (res.empty() && geogSrcIn && vertDstIn &&
+            geogSrcIn->coordinateSystem()->axisList().size() == 3 &&
+            geogSrcIn->datum()) {
+            const auto &authFactory = context.context->getAuthorityFactory();
+            const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
+                authFactory, geogSrcIn->datum().get()));
+            for (const auto &candidate : candidatesSrcGeod) {
+                auto geogCandidate =
+                    util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
+                        candidate);
+                if (geogCandidate &&
+                    geogCandidate->coordinateSystem()->axisList().size() == 2) {
+                    res = findOpsInRegistryDirect(NN_NO_CHECK(geogCandidate),
+                                                  targetCRSIn, context.context);
+                    if (res.empty()) {
+                        res = applyInverse(findOpsInRegistryDirect(
+                            targetCRSIn, NN_NO_CHECK(geogCandidate),
+                            context.context));
+                    }
+                    break;
+                }
+            }
+            return true;
+        }
+        return false;
+    };
+
+    if (geog3DToVertTryThroughGeog2D(geogSrc, vertDst, targetCRS)) {
+        // do nothing
+    } else if (geog3DToVertTryThroughGeog2D(geogDst, vertSrc, sourceCRS)) {
+        res = applyInverse(res);
+    }
+
+    // There's no direct transformation from NAVD88 height to WGS84,
+    // so try to research all transformations from NAVD88 to another
+    // intermediate GeographicCRS.
+    if (res.empty() &&
+        !context.inCreateOperationsGeogToVertWithAlternativeGeog && geogSrc &&
+        vertDst) {
+        res = createOperationsGeogToVertWithAlternativeGeog(sourceCRS,
+                                                            targetCRS, context);
+    } else if (res.empty() &&
+               !context.inCreateOperationsGeogToVertWithAlternativeGeog &&
+               geogDst && vertSrc) {
+        res = applyInverse(createOperationsGeogToVertWithAlternativeGeog(
+            targetCRS, sourceCRS, context));
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
     Private::Context &context, const crs::GeodeticCRS *geodSrc,
@@ -13346,6 +13474,8 @@ getOps(const CoordinateOperationNNPtr &op) {
     return {op};
 }
 
+// ---------------------------------------------------------------------------
+
 static bool useDifferentTransformationsForSameSourceTarget(
     const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) {
     auto subOpsA = getOps(opA);
@@ -13393,6 +13523,63 @@ static bool useDifferentTransformationsForSameSourceTarget(
     return false;
 }
 
+// ---------------------------------------------------------------------------
+
+static crs::GeographicCRSPtr
+getInterpolationGeogCRS(const CoordinateOperationNNPtr &verticalTransform,
+                        const io::DatabaseContextPtr &dbContext) {
+    crs::GeographicCRSPtr interpolationGeogCRS;
+    auto transformationVerticalTransform =
+        dynamic_cast<const Transformation *>(verticalTransform.get());
+    if (transformationVerticalTransform == nullptr) {
+        const auto concat = dynamic_cast<const ConcatenatedOperation *>(
+            verticalTransform.get());
+        if (concat) {
+            const auto &steps = concat->operations();
+            // Is this change of unit + transformation ?
+            if (steps.size() == 2 &&
+                dynamic_cast<const Conversion *>(steps[0].get())) {
+                transformationVerticalTransform =
+                    dynamic_cast<const Transformation *>(steps[1].get());
+            }
+        }
+    }
+    if (transformationVerticalTransform &&
+        !transformationVerticalTransform->hasBallparkTransformation()) {
+        auto interpTransformCRS =
+            transformationVerticalTransform->interpolationCRS();
+        if (interpTransformCRS) {
+            interpolationGeogCRS =
+                std::dynamic_pointer_cast<crs::GeographicCRS>(
+                    interpTransformCRS);
+        } else {
+            // If no explicit interpolation CRS, then
+            // this will be the geographic CRS of the
+            // vertical to geog transformation
+            interpolationGeogCRS =
+                std::dynamic_pointer_cast<crs::GeographicCRS>(
+                    transformationVerticalTransform->targetCRS().as_nullable());
+        }
+    }
+
+    if (interpolationGeogCRS) {
+        if (interpolationGeogCRS->coordinateSystem()->axisList().size() == 3) {
+            // We need to force the interpolation CRS, which
+            // will
+            // frequently be 3D, to 2D to avoid transformations
+            // between source CRS and interpolation CRS to have
+            // 3D terms.
+            interpolationGeogCRS =
+                interpolationGeogCRS->demoteTo2D(std::string(), dbContext)
+                    .as_nullable();
+        }
+    }
+
+    return interpolationGeogCRS;
+}
+
+// ---------------------------------------------------------------------------
+
 void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
     Private::Context &context, const crs::CompoundCRS *compoundSrc,
@@ -13440,7 +13627,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
                 CoordinateOperationContext::GridAvailabilityUse::
                     IGNORE_GRID_AVAILABILITY;
             for (const auto &op : verticalTransforms) {
-                if (!op->identifiers().empty() && dbContext) {
+                if (hasIdentifiers(op) && dbContext) {
                     bool missingGrid = false;
                     if (!ignoreMissingGrids) {
                         const auto gridsNeeded = op->gridsNeeded(dbContext);
@@ -13469,7 +13656,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
                 bool foundRegisteredTransform = false;
                 foundRegisteredTransformWithAllGridsAvailable = false;
                 for (const auto &op : verticalTransformsTmp) {
-                    if (!op->identifiers().empty() && dbContext) {
+                    if (hasIdentifiers(op) && dbContext) {
                         bool missingGrid = false;
                         if (!ignoreMissingGrids) {
                             const auto gridsNeeded = op->gridsNeeded(dbContext);
@@ -13510,43 +13697,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
             cacheHorizToInterpAndInterpToTarget;
 
         for (const auto &verticalTransform : verticalTransforms) {
-            crs::GeographicCRSPtr interpolationGeogCRS;
-            auto transformationVerticalTransform =
-                dynamic_cast<const Transformation *>(verticalTransform.get());
-            if (transformationVerticalTransform &&
-                !transformationVerticalTransform->hasBallparkTransformation()) {
-                auto interpTransformCRS =
-                    transformationVerticalTransform->interpolationCRS();
-                if (interpTransformCRS) {
-                    interpolationGeogCRS =
-                        std::dynamic_pointer_cast<crs::GeographicCRS>(
-                            interpTransformCRS);
-                } else {
-                    // If no explicit interpolation CRS, then
-                    // this will be the geographic CRS of the
-                    // vertical to geog transformation
-                    interpolationGeogCRS =
-                        std::dynamic_pointer_cast<crs::GeographicCRS>(
-                            transformationVerticalTransform->targetCRS()
-                                .as_nullable());
-                }
-            }
-
+            crs::GeographicCRSPtr interpolationGeogCRS =
+                getInterpolationGeogCRS(verticalTransform, dbContext);
             if (interpolationGeogCRS) {
-                if (interpolationGeogCRS->coordinateSystem()
-                        ->axisList()
-                        .size() == 3) {
-                    // We need to force the interpolation CRS, which
-                    // will
-                    // frequently be 3D, to 2D to avoid transformations
-                    // between source CRS and interpolation CRS to have
-                    // 3D terms.
-                    interpolationGeogCRS =
-                        interpolationGeogCRS
-                            ->demoteTo2D(std::string(), dbContext)
-                            .as_nullable();
-                }
-
                 std::vector<CoordinateOperationNNPtr> srcToInterpOps;
                 std::vector<CoordinateOperationNNPtr> interpToTargetOps;
 
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 993b4eeb..d3ea8f26 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -4603,6 +4603,30 @@ std::list<crs::GeodeticCRSNNPtr> AuthorityFactory::createGeodeticCRSFromDatum(
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+std::list<crs::VerticalCRSNNPtr> AuthorityFactory::createVerticalCRSFromDatum(
+    const std::string &datum_auth_name, const std::string &datum_code) const {
+    std::string sql(
+        "SELECT auth_name, code FROM vertical_crs WHERE "
+        "datum_auth_name = ? AND datum_code = ? AND deprecated = 0");
+    ListOfParams params{datum_auth_name, datum_code};
+    if (d->hasAuthorityRestriction()) {
+        sql += " AND auth_name = ?";
+        params.emplace_back(d->authority());
+    }
+    auto sqlRes = d->run(sql, params);
+    std::list<crs::VerticalCRSNNPtr> res;
+    for (const auto &row : sqlRes) {
+        const auto &auth_name = row[0];
+        const auto &code = row[1];
+        res.emplace_back(d->createFactory(auth_name)->createVerticalCRS(code));
+    }
+    return res;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 std::list<crs::GeodeticCRSNNPtr>
 AuthorityFactory::createGeodeticCRSFromEllipsoid(
     const std::string &ellipsoid_auth_name, const std::string &ellipsoid_code,
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 79541d88..7fbbf70e 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -7165,6 +7165,127 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_geogCRS_with_vertical_unit_change) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    // NAD83(2011) + NAVD88 height (ftUS)
+    auto srcObj = createFromUserInput("EPSG:6318+6360",
+                                      authFactory->databaseContext(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+    ASSERT_TRUE(src != nullptr);
+    auto nnSrc = NN_NO_CHECK(src);
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("6319"); // NAD83(2011) 3D
+
+    auto listCompoundToGeog =
+        CoordinateOperationFactory::create()->createOperations(nnSrc, dst,
+                                                               ctxt);
+    ASSERT_TRUE(!listCompoundToGeog.empty());
+
+    // NAD83(2011) + NAVD88 height
+    auto srcObjCompoundVMetre = createFromUserInput(
+        "EPSG:6318+5703", authFactory->databaseContext(), false);
+    auto srcCompoundVMetre = nn_dynamic_pointer_cast<CRS>(srcObjCompoundVMetre);
+    ASSERT_TRUE(srcCompoundVMetre != nullptr);
+    auto listCompoundMetreToGeog =
+        CoordinateOperationFactory::create()->createOperations(
+            NN_NO_CHECK(srcCompoundVMetre), dst, ctxt);
+
+    // Check that we get the same and similar results whether we start from
+    // regular NAVD88 height or its ftUs variant
+    ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size());
+
+    EXPECT_EQ(listCompoundToGeog[0]->nameStr(),
+              "Transformation from NAVD88 height (ftUS) to NAVD88 height + " +
+                  listCompoundMetreToGeog[0]->nameStr());
+    EXPECT_EQ(
+        listCompoundToGeog[0]->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
+                                        authFactory->databaseContext())
+                .get()),
+        replaceAll(listCompoundMetreToGeog[0]->exportToPROJString(
+                       PROJStringFormatter::create(
+                           PROJStringFormatter::Convention::PROJ_5,
+                           authFactory->databaseContext())
+                           .get()),
+                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad",
+                   "+step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad "
+                   "+z_out=m"));
+
+    // Check reverse path
+    auto listGeogToCompound =
+        CoordinateOperationFactory::create()->createOperations(dst, nnSrc,
+                                                               ctxt);
+    EXPECT_EQ(listGeogToCompound.size(), listCompoundToGeog.size());
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    compoundCRS_to_geogCRS_with_vertical_unit_change_and_complex_horizontal_change) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    // NAD83(2011) + NAVD88 height (ftUS)
+    auto srcObj = createFromUserInput("EPSG:6318+6360",
+                                      authFactory->databaseContext(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+    ASSERT_TRUE(src != nullptr);
+    auto nnSrc = NN_NO_CHECK(src);
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("7665"); // WGS84(G1762) 3D
+
+    auto listCompoundToGeog =
+        CoordinateOperationFactory::create()->createOperations(nnSrc, dst,
+                                                               ctxt);
+
+    // NAD83(2011) + NAVD88 height
+    auto srcObjCompoundVMetre = createFromUserInput(
+        "EPSG:6318+5703", authFactory->databaseContext(), false);
+    auto srcCompoundVMetre = nn_dynamic_pointer_cast<CRS>(srcObjCompoundVMetre);
+    ASSERT_TRUE(srcCompoundVMetre != nullptr);
+    auto listCompoundMetreToGeog =
+        CoordinateOperationFactory::create()->createOperations(
+            NN_NO_CHECK(srcCompoundVMetre), dst, ctxt);
+
+    // Check that we get the same and similar results whether we start from
+    // regular NAVD88 height or its ftUs variant
+    ASSERT_EQ(listCompoundToGeog.size(), listCompoundMetreToGeog.size());
+
+    ASSERT_GE(listCompoundToGeog.size(), 1U);
+
+    EXPECT_EQ(listCompoundToGeog[0]->nameStr(),
+              "Transformation from NAVD88 height (ftUS) to NAVD88 height + " +
+                  listCompoundMetreToGeog[0]->nameStr());
+    EXPECT_EQ(
+        listCompoundToGeog[0]->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
+                                        authFactory->databaseContext())
+                .get()),
+        replaceAll(listCompoundMetreToGeog[0]->exportToPROJString(
+                       PROJStringFormatter::create(
+                           PROJStringFormatter::Convention::PROJ_5,
+                           authFactory->databaseContext())
+                           .get()),
+                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad",
+                   "+step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad "
+                   "+z_out=m"));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");