Improve vertical transformation support

When we had a transformation between a compoundCRS and a target geographicCRS,
we did not take into account that in the vertical->other_geog_CRS transformation
we used, the other_geog_CRS was an implicit interpolation CRS. Thus before
doing vertical adjustment, we must go to this interpolation CRS.
The workflow is thus:
source CRS -> interpolation CRS + vertical adjustment + interplation CRS -> target CRS

5 files changed, 443 insertions, 119 deletions

diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 65459058..c07904c8 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -345,6 +345,10 @@ class PROJ_GCC_DLL GeographicCRS : public GeodeticCRS {
            const datum::DatumEnsemblePtr &datumEnsemble,
            const cs::EllipsoidalCSNNPtr &cs);
 
+    PROJ_DLL GeographicCRSNNPtr
+    demoteTo2D(const std::string &newName,
+               const io::DatabaseContextPtr &dbContext) const;
+
     PROJ_DLL static const GeographicCRSNNPtr EPSG_4267; // NAD27
     PROJ_DLL static const GeographicCRSNNPtr EPSG_4269; // NAD83
     PROJ_DLL static const GeographicCRSNNPtr EPSG_4326; // WGS 84 2D
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index b70d3b9d..c9f44ead 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -151,6 +151,7 @@ osgeo::proj::crs::GeodeticCRS::velocityModel() const
 osgeo::proj::crs::GeographicCRS::coordinateSystem() const
 osgeo::proj::crs::GeographicCRS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::datum::GeodeticReferenceFrame> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::EllipsoidalCS> > const&)
 osgeo::proj::crs::GeographicCRS::create(osgeo::proj::util::PropertyMap const&, std::shared_ptr<osgeo::proj::datum::GeodeticReferenceFrame> const&, std::shared_ptr<osgeo::proj::datum::DatumEnsemble> const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::EllipsoidalCS> > const&)
+osgeo::proj::crs::GeographicCRS::demoteTo2D(std::string const&, std::shared_ptr<osgeo::proj::io::DatabaseContext> const&) const
 osgeo::proj::crs::GeographicCRS::~GeographicCRS()
 osgeo::proj::crs::GeographicCRS::is2DPartOf3D(dropbox::oxygen::nn<osgeo::proj::crs::GeographicCRS const*>) const
 osgeo::proj::crs::ParametricCRS::coordinateSystem() const
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index bf3b5146..d4f2903f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -58,7 +58,9 @@
 
 // #define DEBUG
 // #define DEBUG_SORT
-#if defined(DEBUG) || defined(DEBUG_SORT)
+// #define DEBUG_CONCATENATED_OPERATION
+#if defined(DEBUG) || defined(DEBUG_SORT) ||                                   \
+    defined(DEBUG_CONCATENATED_OPERATION)
 #include <iostream>
 #endif
 
@@ -121,6 +123,8 @@ static const char *BALLPARK_VERTICAL_TRANSFORMATION =
 static const char *BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT =
     " (ballpark vertical transformation, without ellipsoid height to vertical "
     "height correction)";
+
+static const std::string AXIS_ORDER_CHANGE_2D_NAME = "axis order change (2D)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -400,7 +404,7 @@ static double getAccuracy(const CoordinateOperationNNPtr &op) {
 
 // ---------------------------------------------------------------------------
 
-// Returns the accuracy of a set of concantenated operations, or -1 if unknown
+// Returns the accuracy of a set of concatenated operations, or -1 if unknown
 static double getAccuracy(const std::vector<CoordinateOperationNNPtr> &ops) {
     double accuracy = -1.0;
     for (const auto &subop : ops) {
@@ -4631,7 +4635,7 @@ ConversionNNPtr Conversion::createAxisOrderReversal(bool is3D) {
                           EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_3D),
                       {}, {});
     } else {
-        return create(createMapNameEPSGCode("axis order change (2D)", 15498),
+        return create(createMapNameEPSGCode(AXIS_ORDER_CHANGE_2D_NAME, 15498),
                       createMethodMapNameEPSGCode(
                           EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_2D),
                       {}, {});
@@ -7421,6 +7425,8 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
                 }
                 if (starts_with(tokens[i], INVERSE_OF)) {
                     name += tokens[i].substr(INVERSE_OF.size());
+                } else if (tokens[i] == AXIS_ORDER_CHANGE_2D_NAME) {
+                    name += tokens[i];
                 } else {
                     name += INVERSE_OF + tokens[i];
                 }
@@ -9405,6 +9411,10 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
         if (i >= 1) {
             if (!compareStepCRS(l_sourceCRS.get(), lastTargetCRS.get())) {
 #ifdef DEBUG_CONCATENATED_OPERATION
+                std::cerr << "Step " << i - 1 << ": "
+                          << operationsIn[i - 1]->nameStr() << std::endl;
+                std::cerr << "Step " << i << ": " << operationsIn[i]->nameStr()
+                          << std::endl;
                 {
                     auto f(io::WKTFormatter::create(
                         io::WKTFormatter::Convention::WKT2_2019));
@@ -9560,7 +9570,7 @@ void ConcatenatedOperation::fixStepsDirection(
             !compareStepCRS(l_sourceCRS.get(), concatOpSourceCRS.get())) {
             throw InvalidOperation("The source CRS of the first step of "
                                    "concatenated operation is not the same "
-                                   "as the source CRS of the concantenated "
+                                   "as the source CRS of the concatenated "
                                    "operation itself");
         }
 
@@ -9569,7 +9579,7 @@ void ConcatenatedOperation::fixStepsDirection(
             !compareStepCRS(l_targetCRS.get(), concatOpTargetCRS.get())) {
             throw InvalidOperation("The target CRS of the last step of "
                                    "concatenated operation is not the same "
-                                   "as the target CRS of the concantenated "
+                                   "as the target CRS of the concatenated "
                                    "operation itself");
         }
     }
@@ -9648,7 +9658,7 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
     auto extent = getExtent(flattenOps, false, emptyIntersection);
     if (checkExtent && emptyIntersection) {
         std::string msg(
-            "empty intersection of area of validity of concantenated "
+            "empty intersection of area of validity of concatenated "
             "operations");
         throw InvalidOperationEmptyIntersection(msg);
     }
@@ -11265,12 +11275,9 @@ createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
         geogSrc->datum()->_isEquivalentTo(
             geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT);
 
-    std::string name(isSameDatum ? NULL_GEOGRAPHIC_OFFSET
-                                 : BALLPARK_GEOGRAPHIC_OFFSET);
-    name += " from ";
-    name += sourceCRS->nameStr();
-    name += " to ";
-    name += targetCRS->nameStr();
+    auto name = buildOpName(isSameDatum ? NULL_GEOGRAPHIC_OFFSET
+                                        : BALLPARK_GEOGRAPHIC_OFFSET,
+                            sourceCRS, targetCRS);
 
     const auto &sourceCRSExtent = getExtent(sourceCRS);
     const auto &targetCRSExtent = getExtent(targetCRS);
@@ -11480,18 +11487,8 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
     auto exportable = util::nn_make_shared<MyPROJStringExportableHorizVertical>(
         horizTransform, verticalTransform, geogDst);
 
-    bool horizTransformIsNoOp = horizTransform->sourceCRS()->_isEquivalentTo(
-        horizTransform->targetCRS().get());
-    if (!horizTransformIsNoOp) {
-        const crs::GeographicCRS *geogSrc =
-            dynamic_cast<const crs::GeographicCRS *>(
-                horizTransform->sourceCRS().get());
-        if (geogSrc) {
-            horizTransformIsNoOp =
-                geogSrc->is2DPartOf3D(NN_NO_CHECK(geogDst.get()));
-        }
-    }
-
+    const bool horizTransformIsNoOp =
+        starts_with(horizTransform->nameStr(), NULL_GEOGRAPHIC_OFFSET);
     if (horizTransformIsNoOp) {
         auto properties = util::PropertyMap();
         properties.set(common::IdentifiedObject::NAME_KEY,
@@ -11541,26 +11538,34 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
     const operation::CoordinateOperationNNPtr &opSrcCRSToGeogCRS,
     const operation::CoordinateOperationNNPtr &verticalTransform,
     const operation::CoordinateOperationNNPtr &opGeogCRStoDstCRS,
-    const crs::GeographicCRSPtr &interpolationGeogCRS) {
+    const crs::GeographicCRSPtr &interpolationGeogCRS, bool checkExtent) {
 
     auto exportable =
         util::nn_make_shared<MyPROJStringExportableHorizVerticalHorizPROJBased>(
             opSrcCRSToGeogCRS, verticalTransform, opGeogCRStoDstCRS,
             interpolationGeogCRS);
 
-    bool dummy = false;
-    auto ops = opSrcCRSToGeogCRS->sourceCRS()->_isEquivalentTo(
-                   opSrcCRSToGeogCRS->targetCRS().get())
-                   ? std::vector<CoordinateOperationNNPtr>{verticalTransform,
-                                                           opGeogCRStoDstCRS}
-                   : std::vector<CoordinateOperationNNPtr>{opSrcCRSToGeogCRS,
-                                                           verticalTransform,
-                                                           opGeogCRStoDstCRS};
+    std::vector<CoordinateOperationNNPtr> ops;
+    if (!starts_with(opSrcCRSToGeogCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+        ops.emplace_back(opSrcCRSToGeogCRS);
+    }
+    ops.emplace_back(verticalTransform);
+    if (!starts_with(opGeogCRStoDstCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+        ops.emplace_back(opGeogCRStoDstCRS);
+    }
+
     bool hasBallparkTransformation = false;
     for (const auto &op : ops) {
         hasBallparkTransformation |= op->hasBallparkTransformation();
     }
-    auto extent = getExtent(ops, true, dummy);
+    bool emptyIntersection = false;
+    auto extent = getExtent(ops, false, emptyIntersection);
+    if (checkExtent && emptyIntersection) {
+        std::string msg(
+            "empty intersection of area of validity of concatenated "
+            "operations");
+        throw InvalidOperationEmptyIntersection(msg);
+    }
     auto properties = util::PropertyMap();
     properties.set(common::IdentifiedObject::NAME_KEY,
                    computeConcatenatedName(ops));
@@ -12311,6 +12316,67 @@ bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult(
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+
+static std::vector<CoordinateOperationNNPtr>
+getOps(const CoordinateOperationNNPtr &op) {
+    auto concatenated = dynamic_cast<const ConcatenatedOperation *>(op.get());
+    if (concatenated)
+        return concatenated->operations();
+    return {op};
+}
+
+static bool useDifferentTransformationsForSameSourceTarget(
+    const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) {
+    auto subOpsA = getOps(opA);
+    auto subOpsB = getOps(opB);
+    for (const auto &subOpA : subOpsA) {
+        if (!dynamic_cast<const Transformation *>(subOpA.get()))
+            continue;
+        if (subOpA->sourceCRS()->nameStr() == "unknown" ||
+            subOpA->targetCRS()->nameStr() == "unknown")
+            continue;
+        for (const auto &subOpB : subOpsB) {
+            if (!dynamic_cast<const Transformation *>(subOpB.get()))
+                continue;
+            if (subOpB->sourceCRS()->nameStr() == "unknown" ||
+                subOpB->targetCRS()->nameStr() == "unknown")
+                continue;
+
+            if (subOpA->sourceCRS()->nameStr() ==
+                    subOpB->sourceCRS()->nameStr() &&
+                subOpA->targetCRS()->nameStr() ==
+                    subOpB->targetCRS()->nameStr()) {
+                if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+                    starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+                    continue;
+                }
+
+                if (!subOpA->isEquivalentTo(subOpB.get())) {
+                    return true;
+                }
+            } else if (subOpA->sourceCRS()->nameStr() ==
+                           subOpB->targetCRS()->nameStr() &&
+                       subOpA->targetCRS()->nameStr() ==
+                           subOpB->sourceCRS()->nameStr()) {
+                if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
+                    starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
+                    continue;
+                }
+
+                if (!subOpA->isEquivalentTo(subOpB->inverse().get())) {
+                    return true;
+                }
+            }
+        }
+    }
+    return false;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 std::vector<CoordinateOperationNNPtr>
 CoordinateOperationFactory::Private::createOperations(
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
@@ -13129,71 +13195,132 @@ CoordinateOperationFactory::Private::createOperations(
                     }
                 }
             }
-            if (!horizTransforms.empty() && !verticalTransforms.empty()) {
-                for (const auto &horizTransform : horizTransforms) {
-                    for (const auto &verticalTransform : verticalTransforms) {
-
-                        crs::GeographicCRSPtr interpolationGeogCRS;
-                        auto transformationVerticalTransform =
-                            dynamic_cast<const Transformation *>(
-                                verticalTransform.get());
-                        if (transformationVerticalTransform) {
-                            auto interpTransformCRS =
-                                transformationVerticalTransform
-                                    ->interpolationCRS();
-                            if (interpTransformCRS) {
-                                auto nn_interpTransformCRS =
-                                    NN_NO_CHECK(interpTransformCRS);
-                                if (dynamic_cast<const crs::GeographicCRS *>(
-                                        nn_interpTransformCRS.get())) {
-                                    interpolationGeogCRS =
-                                        util::nn_dynamic_pointer_cast<
-                                            crs::GeographicCRS>(
-                                            nn_interpTransformCRS);
-                                }
-                            }
+
+            if (horizTransforms.empty() || verticalTransforms.empty()) {
+                return horizTransforms;
+            }
+
+            typedef std::pair<std::vector<CoordinateOperationNNPtr>,
+                              std::vector<CoordinateOperationNNPtr>>
+                PairOfTransforms;
+            std::map<std::string, PairOfTransforms>
+                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());
+                    }
+                }
+
+                if (interpolationGeogCRS) {
+                    if (interpolationGeogCRS->coordinateSystem()
+                            ->axisList()
+                            .size() == 3) {
+                        io::DatabaseContextPtr dbContext;
+                        if (authFactory) {
+                            dbContext = authFactory->databaseContext();
                         }
-                        bool done = false;
-                        if (interpolationGeogCRS &&
-                            (interpolationGeogCRS->_isEquivalentTo(
-                                 srcGeogCRS.get(),
-                                 util::IComparable::Criterion::EQUIVALENT) ||
-                             interpolationGeogCRS->is2DPartOf3D(
-                                 NN_NO_CHECK(srcGeogCRS.get())))) {
-                            auto srcToInterp = createOperations(
+                        // 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;
+
+                    std::string key;
+                    const auto &ids = interpolationGeogCRS->identifiers();
+                    if (!ids.empty()) {
+                        key = (*ids.front()->codeSpace()) + ':' +
+                              ids.front()->code();
+                    }
+                    if (!key.empty()) {
+                        auto iter =
+                            cacheHorizToInterpAndInterpToTarget.find(key);
+                        if (iter == cacheHorizToInterpAndInterpToTarget.end()) {
+                            srcToInterpOps = createOperations(
                                 componentsSrc[0],
                                 NN_NO_CHECK(interpolationGeogCRS), context);
-                            auto interpToCompoundHoriz = createOperations(
-                                NN_NO_CHECK(interpolationGeogCRS),
-                                componentsSrc[0], context);
-                            if (!srcToInterp.empty() &&
-                                !interpToCompoundHoriz.empty()) {
-                                auto op = createHorizVerticalHorizPROJBased(
-                                    sourceCRS, componentsSrc[0],
-                                    srcToInterp.front(), verticalTransform,
-                                    interpToCompoundHoriz.front(),
-                                    interpolationGeogCRS);
-                                done = true;
-                                res.emplace_back(
-                                    ConcatenatedOperation::
-                                        createComputeMetadata(
-                                            {op, horizTransform},
-                                            !allowEmptyIntersection));
-                            }
+                            interpToTargetOps = createOperations(
+                                NN_NO_CHECK(interpolationGeogCRS), targetCRS,
+                                context);
+                            cacheHorizToInterpAndInterpToTarget[key] =
+                                PairOfTransforms(srcToInterpOps,
+                                                 interpToTargetOps);
+                        } else {
+                            srcToInterpOps = iter->second.first;
+                            interpToTargetOps = iter->second.second;
                         }
-                        if (!done) {
-                            auto op = createHorizVerticalPROJBased(
-                                sourceCRS, targetCRS, horizTransform,
-                                verticalTransform);
+                    } else {
+                        srcToInterpOps = createOperations(
+                            componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS),
+                            context);
+                        interpToTargetOps =
+                            createOperations(NN_NO_CHECK(interpolationGeogCRS),
+                                             targetCRS, context);
+                    }
+
+                    for (const auto &srcToInterp : srcToInterpOps) {
+                        for (const auto &interpToTarget : interpToTargetOps) {
+
+                            if (useDifferentTransformationsForSameSourceTarget(
+                                    srcToInterp, interpToTarget)) {
+                                continue;
+                            }
 
-                            res.emplace_back(op);
+                            try {
+                                auto op = createHorizVerticalHorizPROJBased(
+                                    sourceCRS, targetCRS, srcToInterp,
+                                    verticalTransform, interpToTarget,
+                                    interpolationGeogCRS, true);
+                                res.emplace_back(op);
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                                continue;
+                            }
                         }
                     }
+                } else {
+                    // This case is probably only correct if
+                    // verticalTransform and horizTransform are independant
+                    // and in particular that verticalTransform does not
+                    // involve a grid, because of the rather arbitrary order
+                    // horizontal then vertical applied
+                    for (const auto &horizTransform : horizTransforms) {
+                        auto op = createHorizVerticalPROJBased(
+                            sourceCRS, targetCRS, horizTransform,
+                            verticalTransform);
+                        res.emplace_back(op);
+                    }
                 }
-                return res;
-            } else {
-                return horizTransforms;
             }
+
+            return res;
         }
     } else if (compoundSrc && geodDst) {
         auto datum = geodDst->datum();
@@ -13291,10 +13418,15 @@ CoordinateOperationFactory::Private::createOperations(
                     for (const auto &opSrc : opSrcCRSToGeogCRS) {
                         for (const auto &opDst : opGeogCRStoDstCRS) {
 
-                            auto op = createHorizVerticalHorizPROJBased(
-                                sourceCRS, targetCRS, opSrc, verticalTransform,
-                                opDst, interpolationGeogCRS);
-                            res.emplace_back(op);
+                            try {
+                                auto op = createHorizVerticalHorizPROJBased(
+                                    sourceCRS, targetCRS, opSrc,
+                                    verticalTransform, opDst,
+                                    interpolationGeogCRS, true);
+                                res.emplace_back(op);
+                            } catch (
+                                const InvalidOperationEmptyIntersection &) {
+                            }
                         }
                     }
                 }
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 1b7078c3..e4e05166 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -2111,6 +2111,58 @@ GeographicCRSNNPtr GeographicCRS::createEPSG_4807() {
 
 // ---------------------------------------------------------------------------
 
+/** \brief Return a variant of this CRS "demoted" to a 2D one, if not already
+ * the case.
+ *
+ *
+ * @param newName Name of the new CRS. If empty, nameStr() will be used.
+ * @param dbContext Database context to look for potentially already registered
+ *                  2D CRS. May be nullptr.
+ * @return a new CRS demoted to 2D, or the current one if already 2D or not
+ * applicable.
+ * @since 7.0
+ */
+GeographicCRSNNPtr
+GeographicCRS::demoteTo2D(const std::string &newName,
+                          const io::DatabaseContextPtr &dbContext) const {
+
+    const auto &axisList = coordinateSystem()->axisList();
+    if (axisList.size() == 3) {
+        const auto &l_identifiers = identifiers();
+        // First check if there is a Geographic 3D CRS in the database
+        // of the same name.
+        // This is the common practice in the EPSG dataset.
+        if (dbContext && l_identifiers.size() == 1) {
+            auto authFactory = io::AuthorityFactory::create(
+                NN_NO_CHECK(dbContext), *(l_identifiers[0]->codeSpace()));
+            auto res = authFactory->createObjectsFromName(
+                nameStr(),
+                {io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS}, false);
+            if (!res.empty()) {
+                const auto &firstRes = res.front();
+                auto firstResAsGeogCRS =
+                    util::nn_dynamic_pointer_cast<GeographicCRS>(firstRes);
+                if (firstResAsGeogCRS &&
+                    firstResAsGeogCRS->is2DPartOf3D(NN_NO_CHECK(this))) {
+                    return NN_NO_CHECK(firstResAsGeogCRS);
+                }
+            }
+        }
+
+        auto cs = cs::EllipsoidalCS::create(util::PropertyMap(), axisList[0],
+                                            axisList[1]);
+        return GeographicCRS::create(
+            util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                    !newName.empty() ? newName : nameStr()),
+            datum(), datumEnsemble(), cs);
+    }
+
+    return NN_NO_CHECK(std::dynamic_pointer_cast<GeographicCRS>(
+        shared_from_this().as_nullable()));
+}
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 void GeographicCRS::addAngularUnitConvertAndAxisSwap(
     io::PROJStringFormatter *formatter) const {
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index c3ed371e..ddbb1cc9 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -6115,10 +6115,10 @@ TEST(operation, WGS84_G1762_to_compoundCRS_with_bound_vertCRS) {
         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 WGS 84 to WGS 84 (G1762) + "
+              "Inverse of unknown to WGS84 ellipsoidal height + "
               "Inverse of NAD83 to WGS 84 (1) + "
-              "Inverse of axis order change (2D)");
+              "axis order change (2D)");
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline "
               "+step +proj=axisswap +order=2,1 "
@@ -6465,15 +6465,19 @@ TEST(operation, compoundCRS_with_boundGeogCRS_and_boundVerticalCRS_to_geogCRS) {
     // Not completely sure the order of horizontal and vertical operations
     // makes sense
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
-              "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 "
-              "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
-              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
-              "+inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
-              "+proj=vgridshift +grids=egm08_25.gtx +multiplier=1 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1");
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris "
+              "+step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign "
+              "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+              "+convention=position_vector "
+              "+step +inv +proj=cart +ellps=WGS84 "
+              "+step +proj=pop +v_3 "
+              "+step +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 
     auto grids = op->gridsNeeded(DatabaseContext::create());
     EXPECT_EQ(grids.size(), 1U);
@@ -6503,13 +6507,17 @@ TEST(operation, compoundCRS_with_boundProjCRS_and_boundVerticalCRS_to_geogCRS) {
     // Not completely sure the order of horizontal and vertical operations
     // makes sense
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
-              "+pm=paris +step +proj=push +v_3 +step +proj=cart "
-              "+ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
-              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-              "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=vgridshift "
-              "+grids=egm08_25.gtx +multiplier=1 +step +proj=unitconvert "
-              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=31 +ellps=clrk80ign +pm=paris "
+              "+step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign "
+              "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+              "+convention=position_vector "
+              "+step +inv +proj=cart +ellps=WGS84 "
+              "+step +proj=pop +v_3 "
+              "+step +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 
     auto opInverse = CoordinateOperationFactory::create()->createOperation(
         GeographicCRS::EPSG_4979, compound);
@@ -6890,8 +6898,9 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
                   "+proj=axisswap +order=2,1");
     }
 
-    // CompoundCRS to Geog3DCRS, with same vertical unit, but without
-    // ellipsoid height <--> vertical height correction
+    // CompoundCRS to Geog3DCRS, with same vertical unit, and with
+    // direct ellipsoid height <--> vertical height correction and
+    // direct horizontal transform (no-op here)
     {
         auto ctxt =
             CoordinateOperationContext::create(authFactory, nullptr, 0.0);
@@ -6905,16 +6914,142 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
             ctxt);
         ASSERT_GE(list.size(), 1U);
         EXPECT_EQ(list[0]->nameStr(),
-                  "NAD83(NSRS2007) to WGS 84 (1) + "
-                  "Inverse of NAD83(2011) to NAVD88 height (1)");
+                  "Inverse of NAD83(NSRS2007) to NAVD88 height (1) + "
+                  "NAD83(NSRS2007) to WGS 84 (1)");
         EXPECT_EQ(list[0]->exportToPROJString(
                       PROJStringFormatter::create(
                           PROJStringFormatter::Convention::PROJ_5,
                           authFactory->databaseContext())
                           .get()),
-                  "+proj=pipeline +step +proj=axisswap +order=2,1 "
+                  "+proj=pipeline "
+                  "+step +proj=axisswap +order=2,1 "
+                  "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+                  // Inv here since the grid is not known
+                  "+step +inv +proj=vgridshift +grids=geoid09_conus.bin "
+                  "+multiplier=1 "
+                  "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+                  "+step +proj=axisswap +order=2,1");
+    }
+
+    // CompoundCRS to Geog3DCRS, with same vertical unit, and with
+    // ellipsoid height <--> vertical height correction that requires a
+    // horizontal adjustment before and after (which is empty in practice here
+    // as NAD83 to NAD83(2011) is no-op)
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        ctxt->setSpatialCriterion(
+            CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+        ctxt->setGridAvailabilityUse(
+            CoordinateOperationContext::GridAvailabilityUse::
+                IGNORE_GRID_AVAILABILITY);
+        // NAD83 + NAVD88 height
+        auto srcObj = createFromUserInput(
+            "EPSG:4269+5703", authFactory->databaseContext(), false);
+        auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+        ASSERT_TRUE(src != nullptr);
+        auto nnSrc = NN_NO_CHECK(src);
+
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            nnSrc,
+            authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
+            ctxt);
+        ASSERT_GE(list.size(), 2U);
+
+        EXPECT_EQ(list[0]->nameStr(),
+                  "NAD83 to NAD83(2011) (1) + "
+                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
+                  "Inverse of NAD83 to NAD83(2011) (1) + "
+                  "NAD83 to WGS 84 (1)");
+        EXPECT_EQ(list[0]->exportToPROJString(
+                      PROJStringFormatter::create(
+                          PROJStringFormatter::Convention::PROJ_5,
+                          authFactory->databaseContext())
+                          .get()),
+                  "+proj=pipeline "
+                  "+step +proj=axisswap +order=2,1 "
+                  "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+multiplier=1 "
+                  "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+                  "+step +proj=axisswap +order=2,1");
+
+        // Shows vertical step, and then horizontal step
+        EXPECT_EQ(list[1]->nameStr(),
+                  "NAD83 to NAD83(2011) (1) + "
+                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
+                  "Inverse of NAD83 to NAD83(2011) (1) + "
+                  "NAD83 to WGS 84 (18)");
+        EXPECT_EQ(list[1]->exportToPROJString(
+                      PROJStringFormatter::create(
+                          PROJStringFormatter::Convention::PROJ_5,
+                          authFactory->databaseContext())
+                          .get()),
+                  "+proj=pipeline "
+                  "+step +proj=axisswap +order=2,1 "
+                  "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+multiplier=1 "
+                  "+step +proj=hgridshift +grids=FL "
+                  "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+                  "+step +proj=axisswap +order=2,1");
+    }
+
+    // Another variation, but post horizontal adjustment is in two steps
+    {
+        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
+        auto srcObj = createFromUserInput(
+            "EPSG:6318+5703", authFactory->databaseContext(), false);
+        auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+        ASSERT_TRUE(src != nullptr);
+        auto nnSrc = NN_NO_CHECK(src);
+
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            nnSrc,
+            authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
+            ctxt);
+        ASSERT_GE(list.size(), 2U);
+
+        EXPECT_EQ(list[0]->nameStr(),
+                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
+                  "Inverse of NAD83 to NAD83(2011) (1) + "
+                  "NAD83 to WGS 84 (1)");
+        EXPECT_EQ(list[0]->exportToPROJString(
+                      PROJStringFormatter::create(
+                          PROJStringFormatter::Convention::PROJ_5,
+                          authFactory->databaseContext())
+                          .get()),
+                  "+proj=pipeline "
+                  "+step +proj=axisswap +order=2,1 "
+                  "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+multiplier=1 "
+                  "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+                  "+step +proj=axisswap +order=2,1");
+
+        // Shows vertical step, and then horizontal step
+        EXPECT_EQ(list[1]->nameStr(),
+                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
+                  "Inverse of NAD83 to NAD83(2011) (1) + "
+                  "NAD83 to WGS 84 (18)");
+        EXPECT_EQ(list[1]->exportToPROJString(
+                      PROJStringFormatter::create(
+                          PROJStringFormatter::Convention::PROJ_5,
+                          authFactory->databaseContext())
+                          .get()),
+                  "+proj=pipeline "
+                  "+step +proj=axisswap +order=2,1 "
                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
-                  "+step +proj=vgridshift +grids=g2012bu0.gtx +multiplier=1 "
+                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+multiplier=1 "
+                  "+step +proj=hgridshift +grids=FL "
                   "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
                   "+step +proj=axisswap +order=2,1");
     }
@@ -7002,7 +7137,7 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
             "+ellps=bessel "
             "+step +proj=vgridshift +grids=naptrans2008.gtx +multiplier=1 "
             "+step +proj=hgridshift +grids=rdtrans2008.gsb "
-            "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
             "+step +proj=axisswap +order=2,1");
     }