Merge remote-tracking branch 'osgeo/master'

6 files changed, 1601 insertions, 1070 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index f5f7ba55..337de313 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -498,7 +498,7 @@ template <class T> static PROJ_STRING_LIST to_string_list(T &&set) {
  * proj_string_list_destroy().
  * @param out_grammar_errors Pointer to a PROJ_STRING_LIST object, or NULL.
  * If provided, *out_grammar_errors will contain a list of errors regarding the
- * WKT grammaer. It must be freed with proj_string_list_destroy().
+ * WKT grammar. It must be freed with proj_string_list_destroy().
  * @return Object that must be unreferenced with proj_destroy(), or NULL in
  * case of error.
  */
@@ -522,6 +522,7 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt,
         if (dbContext) {
             parser.attachDatabaseContext(NN_NO_CHECK(dbContext));
         }
+        parser.setStrict(false);
         for (auto iter = options; iter && iter[0]; ++iter) {
             const char *value;
             if ((value = getOptionValue(*iter, "STRICT="))) {
@@ -536,10 +537,19 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt,
         auto obj = nn_dynamic_pointer_cast<IdentifiedObject>(
             parser.createFromWKT(wkt));
 
+        std::vector<std::string> warningsFromParsing;
         if (out_grammar_errors) {
-            auto warnings = parser.warningList();
-            if (!warnings.empty()) {
-                *out_grammar_errors = to_string_list(warnings);
+            auto rawWarnings = parser.warningList();
+            std::vector<std::string> grammarWarnings;
+            for (const auto &msg : rawWarnings) {
+                if (msg.find("Default it to") != std::string::npos) {
+                    warningsFromParsing.push_back(msg);
+                } else {
+                    grammarWarnings.push_back(msg);
+                }
+            }
+            if (!grammarWarnings.empty()) {
+                *out_grammar_errors = to_string_list(grammarWarnings);
             }
         }
 
@@ -548,6 +558,8 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt,
             if (derivedCRS) {
                 auto warnings =
                     derivedCRS->derivingConversionRef()->validateParameters();
+                warnings.insert(warnings.end(), warningsFromParsing.begin(),
+                                warningsFromParsing.end());
                 if (!warnings.empty()) {
                     *out_warnings = to_string_list(warnings);
                 }
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 0f0e216c..9e3306ba 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -3197,7 +3197,7 @@ ConversionNNPtr Conversion::createBonne(const util::PropertyMap &properties,
  * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::9834)
  *
  * \warning The PROJ cea computation code would select the ellipsoidal form if
- * a non-spherical ellipsoid is used for the base GeographicalCRS.
+ * a non-spherical ellipsoid is used for the base GeographicCRS.
  *
  * @param properties See \ref general_properties of the conversion. If the name
  * is not provided, it is automatically set.
@@ -4661,6 +4661,26 @@ Conversion::createChangeVerticalUnit(const util::PropertyMap &properties,
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a conversion based on the Height Depth Reversal
+ * method.
+ *
+ * This method is defined as [EPSG:1068]
+ * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::1068)
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @return a new Conversion.
+ * @since 7.0
+ */
+ConversionNNPtr
+Conversion::createHeightDepthReversal(const util::PropertyMap &properties) {
+    return create(properties, createMethodMapNameEPSGCode(
+                                  EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL),
+                  {}, {});
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instantiate a conversion based on the Axis order reversal method
  *
  * This swaps the longitude, latitude axis.
@@ -4927,6 +4947,14 @@ CoordinateOperationNNPtr Conversion::inverse() const {
         return conv;
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) {
+
+        auto conv = createHeightDepthReversal(
+            createPropertiesForInverse(this, false, false));
+        conv->setCRSs(this, true);
+        return conv;
+    }
+
     return InverseConversion::create(NN_NO_CHECK(
         util::nn_dynamic_pointer_cast<Conversion>(shared_from_this())));
 }
@@ -5808,9 +5836,12 @@ void Conversion::_exportToPROJString(
         methodEPSGCode == EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION;
     const bool isGeographicGeocentric =
         methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC;
+    const bool isHeightDepthReversal =
+        methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL;
     const bool applySourceCRSModifiers =
         !isZUnitConversion && !isAffineParametric &&
-        !isAxisOrderReversal(methodEPSGCode) && !isGeographicGeocentric;
+        !isAxisOrderReversal(methodEPSGCode) && !isGeographicGeocentric &&
+        !isHeightDepthReversal;
     bool applyTargetCRSModifiers = applySourceCRSModifiers;
 
     auto l_sourceCRS = sourceCRS();
@@ -6056,24 +6087,31 @@ void Conversion::_exportToPROJString(
                 if (!param->proj_name) {
                     continue;
                 }
-                auto value =
+                const auto value =
                     parameterValueMeasure(param->wkt2_name, param->epsg_code);
+                double valueConverted = 0;
+                if (value == nullMeasure) {
+                    // Deal with missing values. In an ideal world, this would
+                    // not happen
+                    if (param->epsg_code ==
+                        EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) {
+                        valueConverted = 1.0;
+                    }
+                } else if (param->unit_type ==
+                           common::UnitOfMeasure::Type::ANGULAR) {
+                    valueConverted =
+                        value.convertToUnit(common::UnitOfMeasure::DEGREE);
+                } else {
+                    valueConverted = value.getSIValue();
+                }
+
                 if (mapping->epsg_code ==
                         EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP &&
                     strcmp(param->proj_name, "lat_1") == 0) {
-                    formatter->addParam(
-                        param->proj_name,
-                        value.convertToUnit(common::UnitOfMeasure::DEGREE));
-                    formatter->addParam(
-                        "lat_0",
-                        value.convertToUnit(common::UnitOfMeasure::DEGREE));
-                } else if (param->unit_type ==
-                           common::UnitOfMeasure::Type::ANGULAR) {
-                    formatter->addParam(
-                        param->proj_name,
-                        value.convertToUnit(common::UnitOfMeasure::DEGREE));
+                    formatter->addParam(param->proj_name, valueConverted);
+                    formatter->addParam("lat_0", valueConverted);
                 } else {
-                    formatter->addParam(param->proj_name, value.getSIValue());
+                    formatter->addParam(param->proj_name, valueConverted);
                 }
             }
 
@@ -7876,6 +7914,17 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                 coordinateOperationAccuracies()));
     }
 
+#ifdef notdef
+    // We dont need that currently, but we might...
+    if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) {
+        return d->registerInv(
+            shared_from_this(),
+            createHeightDepthReversal(
+                createPropertiesForInverse(this, false, false), l_targetCRS,
+                l_sourceCRS, coordinateOperationAccuracies()));
+    }
+#endif
+
     return InverseTransformation::create(NN_NO_CHECK(
         util::nn_dynamic_pointer_cast<Transformation>(shared_from_this())));
 }
@@ -9389,6 +9438,12 @@ bool SingleOperation::exportToPROJStringGeneric(
         return true;
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) {
+        formatter->addStep("axisswap");
+        formatter->addParam("order", "1,2,-3");
+        return true;
+    }
+
     return false;
 }
 
@@ -10297,7 +10352,8 @@ struct CoordinateOperationFactory::Private {
         const CoordinateOperationContextNNPtr &context;
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
         bool inCreateOperationsThroughPreferredHub = false;
-        bool inCreateOperationsGeogToVertWithIntermediate = false;
+        bool inCreateOperationsGeogToVertWithAlternativeGeog = false;
+        bool inCreateOperationsGeogToVertWithIntermediateVert = false;
         bool skipHorizontalTransformation = false;
 
         Context(const crs::CRSNNPtr &sourceCRSIn,
@@ -10312,6 +10368,103 @@ struct CoordinateOperationFactory::Private {
                      const crs::CRSNNPtr &targetCRS, Context &context);
 
   private:
+    static constexpr bool allowEmptyIntersection = true;
+
+    static void createOperationsFromProj4Ext(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static bool createOperationsFromDatabase(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::GeodeticCRS *geodSrc,
+        const crs::GeodeticCRS *geodDst, const crs::GeographicCRS *geogSrc,
+        const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
+        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,
+        const crs::GeodeticCRS *geodDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsDerivedTo(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::DerivedCRS *derivedSrc,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsBoundToGeog(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::BoundCRS *boundSrc,
+        const crs::GeographicCRS *geogDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsBoundToVert(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::BoundCRS *boundSrc,
+        const crs::VerticalCRS *vertDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsVertToVert(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::VerticalCRS *vertSrc,
+        const crs::VerticalCRS *vertDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsVertToGeog(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::VerticalCRS *vertSrc,
+        const crs::GeographicCRS *geogDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsBoundToBound(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::BoundCRS *boundSrc,
+        const crs::BoundCRS *boundDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsCompoundToGeog(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::CompoundCRS *compoundSrc,
+        const crs::GeographicCRS *geogDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsCompoundToGeod(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::CompoundCRS *compoundSrc,
+        const crs::GeodeticCRS *geodDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsCompoundToCompound(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::CompoundCRS *compoundSrc,
+        const crs::CompoundCRS *compoundDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
+    static void createOperationsBoundToCompound(
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        Private::Context &context, const crs::BoundCRS *boundSrc,
+        const crs::CompoundCRS *compoundDst,
+        std::vector<CoordinateOperationNNPtr> &res);
+
     static std::vector<CoordinateOperationNNPtr> createOperationsGeogToGeog(
         std::vector<CoordinateOperationNNPtr> &res,
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
@@ -10329,11 +10482,6 @@ struct CoordinateOperationFactory::Private {
         const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst,
         Context &context);
 
-    static std::vector<CoordinateOperationNNPtr>
-    createOperationsGeogToVertWithIntermediate(const crs::CRSNNPtr &sourceCRS,
-                                               const crs::CRSNNPtr &targetCRS,
-                                               Context &context);
-
     static bool
     hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> &res,
                              const Context &context);
@@ -11544,6 +11692,9 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
 
 MyPROJStringExportableHorizVerticalHorizPROJBased::
     ~MyPROJStringExportableHorizVerticalHorizPROJBased() = default;
+
+//! @endcond
+
 } // namespace operation
 NS_PROJ_END
 
@@ -11558,6 +11709,8 @@ template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableHorizVer
 NS_PROJ_START
 namespace operation {
 
+//! @cond Doxygen_Suppress
+
 // ---------------------------------------------------------------------------
 
 static std::string buildTransfName(const std::string &srcName,
@@ -11711,7 +11864,6 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
 
     assert(sourceCRS.get() == geogSrc);
     assert(targetCRS.get() == geogDst);
-    const bool allowEmptyIntersection = true;
 
     const auto &src_pm = geogSrc->primeMeridian()->longitude();
     const auto &dst_pm = geogDst->primeMeridian()->longitude();
@@ -11885,12 +12037,8 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
     return res;
 }
 
-//! @endcond
-
 // ---------------------------------------------------------------------------
 
-//! @cond Doxygen_Suppress
-
 static bool hasIdentifiers(const CoordinateOperationNNPtr &op) {
     if (!op->identifiers().empty()) {
         return true;
@@ -11905,12 +12053,9 @@ static bool hasIdentifiers(const CoordinateOperationNNPtr &op) {
     }
     return false;
 }
-//! @endcond
 
 // ---------------------------------------------------------------------------
 
-//! @cond Doxygen_Suppress
-
 static std::vector<crs::CRSNNPtr>
 findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
                              const datum::GeodeticReferenceFrame *datum) {
@@ -12032,7 +12177,6 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
                objectAsStr(sourceCRS.get()) + "," +
                objectAsStr(targetCRS.get()) + ")");
 #endif
-    const bool allowEmptyIntersection = true;
 
     struct CreateOperationsWithDatumPivotAntiRecursion {
         Context &context;
@@ -12298,7 +12442,6 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub(
         }
     }
 
-    const bool allowEmptyIntersection = true;
     for (const auto &intermCRS : candidatesIntermCRS) {
 #ifdef DEBUG
         EnterDebugLevel loopLevel;
@@ -12327,51 +12470,6 @@ void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub(
 
 // ---------------------------------------------------------------------------
 
-std::vector<CoordinateOperationNNPtr>
-CoordinateOperationFactory::Private::createOperationsGeogToVertWithIntermediate(
-    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.inCreateOperationsGeogToVertWithIntermediate);
-            context.inCreateOperationsGeogToVertWithIntermediate = true;
-        }
-
-        ~AntiRecursionGuard() {
-            context.inCreateOperationsGeogToVertWithIntermediate = 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) {
@@ -12390,12 +12488,8 @@ createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
             sourceCRS, targetCRS, 0.0, 0.0, 0.0, {}));
 }
 
-//! @endcond
-
 // ---------------------------------------------------------------------------
 
-//! @cond Doxygen_Suppress
-
 bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult(
     const std::vector<CoordinateOperationNNPtr> &res, const Context &context) {
     auto resTmp = FilterResults(res, context.context, context.sourceCRS,
@@ -12410,72 +12504,8 @@ bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult(
     return false;
 }
 
-//! @endcond
-
 // ---------------------------------------------------------------------------
 
-//! @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,
@@ -12489,7 +12519,6 @@ CoordinateOperationFactory::Private::createOperations(
 #endif
 
     std::vector<CoordinateOperationNNPtr> res;
-    const bool allowEmptyIntersection = true;
 
     auto boundSrc = dynamic_cast<const crs::BoundCRS *>(sourceCRS.get());
     auto boundDst = dynamic_cast<const crs::BoundCRS *>(targetCRS.get());
@@ -12501,49 +12530,8 @@ CoordinateOperationFactory::Private::createOperations(
                                      ? boundDst->baseCRS()->getExtensionProj4()
                                      : targetCRS->getExtensionProj4();
     if (!sourceProj4Ext.empty() || !targetProj4Ext.empty()) {
-
-        auto sourceProjExportable =
-            dynamic_cast<const io::IPROJStringExportable *>(
-                boundSrc ? boundSrc : sourceCRS.get());
-        auto targetProjExportable =
-            dynamic_cast<const io::IPROJStringExportable *>(
-                boundDst ? boundDst : targetCRS.get());
-        if (!sourceProjExportable) {
-            throw InvalidOperation("Source CRS is not PROJ exportable");
-        }
-        if (!targetProjExportable) {
-            throw InvalidOperation("Target CRS is not PROJ exportable");
-        }
-        auto projFormatter = io::PROJStringFormatter::create();
-        projFormatter->setCRSExport(true);
-        projFormatter->setLegacyCRSToCRSContext(true);
-        projFormatter->startInversion();
-        sourceProjExportable->_exportToPROJString(projFormatter.get());
-        auto geogSrc =
-            dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
-        if (geogSrc) {
-            auto tmpFormatter = io::PROJStringFormatter::create();
-            geogSrc->addAngularUnitConvertAndAxisSwap(tmpFormatter.get());
-            projFormatter->ingestPROJString(tmpFormatter->toString());
-        }
-
-        projFormatter->stopInversion();
-
-        targetProjExportable->_exportToPROJString(projFormatter.get());
-        auto geogDst =
-            dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
-        if (geogDst) {
-            auto tmpFormatter = io::PROJStringFormatter::create();
-            geogDst->addAngularUnitConvertAndAxisSwap(tmpFormatter.get());
-            projFormatter->ingestPROJString(tmpFormatter->toString());
-        }
-
-        const auto PROJString = projFormatter->toString();
-        auto properties = util::PropertyMap().set(
-            common::IdentifiedObject::NAME_KEY,
-            buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
-        res.emplace_back(SingleOperation::createPROJBased(
-            properties, PROJString, sourceCRS, targetCRS, {}));
+        createOperationsFromProj4Ext(sourceCRS, targetCRS, boundSrc, boundDst,
+                                     res);
         return res;
     }
 
@@ -12566,102 +12554,215 @@ CoordinateOperationFactory::Private::createOperations(
          !derivedDst->baseCRS()->_isEquivalentTo(
              sourceCRS.get(), util::IComparable::Criterion::EQUIVALENT))) {
 
-        bool doFilterAndCheckPerfectOp = true;
-        res = findOpsInRegistryDirect(sourceCRS, targetCRS, context.context);
-        if (!sourceCRS->_isEquivalentTo(targetCRS.get())) {
-            auto resFromInverse = applyInverse(
-                findOpsInRegistryDirect(targetCRS, sourceCRS, context.context));
-            res.insert(res.end(), resFromInverse.begin(), resFromInverse.end());
+        if (createOperationsFromDatabase(sourceCRS, targetCRS, context, geodSrc,
+                                         geodDst, geogSrc, geogDst, vertSrc,
+                                         vertDst, res)) {
+            return res;
+        }
+    }
 
-            // If we get at least a result with perfect accuracy, do not
-            // bother generating synthetic transforms.
-            if (hasPerfectAccuracyResult(res, context)) {
-                return res;
-            }
+    // Special case if both CRS are geodetic
+    if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
+        createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
+                                   geodDst, res);
+        return res;
+    }
 
-            doFilterAndCheckPerfectOp = false;
-
-            bool sameGeodeticDatum = false;
-            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.
-            if (res.empty() && geogSrc && vertDst &&
-                geogSrc->coordinateSystem()->axisList().size() == 3 &&
-                geogSrc->datum()) {
-                const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
-                    authFactory, geogSrc->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),
-                                                    targetCRS, context.context);
-                        if (res.empty()) {
-                            res = applyInverse(findOpsInRegistryDirect(
-                                targetCRS, NN_NO_CHECK(geogCandidate),
-                                context.context));
-                        }
-                        break;
-                    }
-                }
-            } else if (res.empty() && geogDst && vertSrc &&
-                       geogDst->coordinateSystem()->axisList().size() == 3 &&
-                       geogDst->datum()) {
-                const auto candidatesDstGeod(findCandidateGeodCRSForDatum(
-                    authFactory, geogDst->datum().get()));
-                for (const auto &candidate : candidatesDstGeod) {
-                    auto geogCandidate =
-                        util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
-                            candidate);
-                    if (geogCandidate &&
-                        geogCandidate->coordinateSystem()->axisList().size() ==
-                            2) {
-                        res = findOpsInRegistryDirect(
-                            sourceCRS, NN_NO_CHECK(geogCandidate),
-                            context.context);
-                        if (res.empty()) {
-                            res = applyInverse(findOpsInRegistryDirect(
-                                NN_NO_CHECK(geogCandidate), sourceCRS,
-                                context.context));
-                        }
-                        break;
-                    }
-                }
-            }
+    // If the source is a derived CRS, then chain the inverse of its
+    // deriving conversion, with transforms from its baseCRS to the
+    // targetCRS
+    if (derivedSrc) {
+        createOperationsDerivedTo(sourceCRS, targetCRS, context, derivedSrc,
+                                  res);
+        return 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.inCreateOperationsGeogToVertWithIntermediate &&
-                geogSrc && vertDst) {
-                res = createOperationsGeogToVertWithIntermediate(
-                    sourceCRS, targetCRS, context);
-            } else if (res.empty() &&
-                       !context.inCreateOperationsGeogToVertWithIntermediate &&
-                       geogDst && vertSrc) {
-                res = applyInverse(createOperationsGeogToVertWithIntermediate(
-                    targetCRS, sourceCRS, context));
-            }
+    // reverse of previous case
+    if (derivedDst) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    if (boundSrc && geogDst) {
+        createOperationsBoundToGeog(sourceCRS, targetCRS, context, boundSrc,
+                                    geogDst, res);
+        return res;
+    }
+
+    // reverse of previous case
+    if (geogSrc && boundDst) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    // vertCRS (as boundCRS with transformation to target vertCRS) to
+    // vertCRS
+    if (boundSrc && vertDst) {
+        createOperationsBoundToVert(sourceCRS, targetCRS, context, boundSrc,
+                                    vertDst, res);
+        return res;
+    }
+
+    // reverse of previous case
+    if (boundDst && vertSrc) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    if (vertSrc && vertDst) {
+        createOperationsVertToVert(sourceCRS, targetCRS, context, vertSrc,
+                                   vertDst, res);
+        return res;
+    }
+
+    // A bit odd case as we are comparing apples to oranges, but in case
+    // the vertical unit differ, do something useful.
+    if (vertSrc && geogDst) {
+        createOperationsVertToGeog(sourceCRS, targetCRS, context, vertSrc,
+                                   geogDst, res);
+        return res;
+    }
+
+    // reverse of previous case
+    if (vertDst && geogSrc) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    // boundCRS to boundCRS using the same geographic hubCRS
+    if (boundSrc && boundDst) {
+        createOperationsBoundToBound(sourceCRS, targetCRS, context, boundSrc,
+                                     boundDst, res);
+        return res;
+    }
+
+    auto compoundSrc = dynamic_cast<crs::CompoundCRS *>(sourceCRS.get());
+    // Order of comparison between the geogDst vs geodDst is impotant
+    if (compoundSrc && geogDst) {
+        createOperationsCompoundToGeog(sourceCRS, targetCRS, context,
+                                       compoundSrc, geogDst, res);
+        return res;
+    } else if (compoundSrc && geodDst) {
+        createOperationsCompoundToGeod(sourceCRS, targetCRS, context,
+                                       compoundSrc, geodDst, res);
+        return res;
+    }
+
+    // reverse of previous cases
+    auto compoundDst = dynamic_cast<const crs::CompoundCRS *>(targetCRS.get());
+    if (geodSrc && compoundDst) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    if (compoundSrc && compoundDst) {
+        createOperationsCompoundToCompound(sourceCRS, targetCRS, context,
+                                           compoundSrc, compoundDst, res);
+        return res;
+    }
+
+    // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to
+    // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx
+    // +type=crs'
+    if (boundSrc && compoundDst) {
+        createOperationsBoundToCompound(sourceCRS, targetCRS, context, boundSrc,
+                                        compoundDst, res);
+        return res;
+    }
+
+    // reverse of previous case
+    if (boundDst && compoundSrc) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
+    return res;
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsFromProj4Ext(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    auto sourceProjExportable = dynamic_cast<const io::IPROJStringExportable *>(
+        boundSrc ? boundSrc : sourceCRS.get());
+    auto targetProjExportable = dynamic_cast<const io::IPROJStringExportable *>(
+        boundDst ? boundDst : targetCRS.get());
+    if (!sourceProjExportable) {
+        throw InvalidOperation("Source CRS is not PROJ exportable");
+    }
+    if (!targetProjExportable) {
+        throw InvalidOperation("Target CRS is not PROJ exportable");
+    }
+    auto projFormatter = io::PROJStringFormatter::create();
+    projFormatter->setCRSExport(true);
+    projFormatter->setLegacyCRSToCRSContext(true);
+    projFormatter->startInversion();
+    sourceProjExportable->_exportToPROJString(projFormatter.get());
+    auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+    if (geogSrc) {
+        auto tmpFormatter = io::PROJStringFormatter::create();
+        geogSrc->addAngularUnitConvertAndAxisSwap(tmpFormatter.get());
+        projFormatter->ingestPROJString(tmpFormatter->toString());
+    }
+
+    projFormatter->stopInversion();
+
+    targetProjExportable->_exportToPROJString(projFormatter.get());
+    auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+    if (geogDst) {
+        auto tmpFormatter = io::PROJStringFormatter::create();
+        geogDst->addAngularUnitConvertAndAxisSwap(tmpFormatter.get());
+        projFormatter->ingestPROJString(tmpFormatter->toString());
+    }
+
+    const auto PROJString = projFormatter->toString();
+    auto properties = util::PropertyMap().set(
+        common::IdentifiedObject::NAME_KEY,
+        buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
+    res.emplace_back(SingleOperation::createPROJBased(
+        properties, PROJString, sourceCRS, targetCRS, {}));
+}
+
+// ---------------------------------------------------------------------------
+
+bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::GeodeticCRS *geodSrc,
+    const crs::GeodeticCRS *geodDst, const crs::GeographicCRS *geogSrc,
+    const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
+    const crs::VerticalCRS *vertDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    bool doFilterAndCheckPerfectOp = true;
+    res = findOpsInRegistryDirect(sourceCRS, targetCRS, context.context);
+    if (!sourceCRS->_isEquivalentTo(targetCRS.get())) {
+        auto resFromInverse = applyInverse(
+            findOpsInRegistryDirect(targetCRS, sourceCRS, context.context));
+        res.insert(res.end(), resFromInverse.begin(), resFromInverse.end());
+
+        // If we get at least a result with perfect accuracy, do not
+        // bother generating synthetic transforms.
+        if (hasPerfectAccuracyResult(res, context)) {
+            return true;
+        }
+
+        doFilterAndCheckPerfectOp = false;
+
+        bool sameGeodeticDatum = false;
+
+        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);
 
             if (res.empty() && !sameGeodeticDatum &&
                 !context.inCreateOperationsWithDatumPivotAntiRecursion &&
-                geodSrc && geodDst) {
+                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
@@ -12671,284 +12772,533 @@ CoordinateOperationFactory::Private::createOperations(
                 // 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) {
-                    createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
-                                                   geodSrc, geodDst, context);
-                    doFilterAndCheckPerfectOp = !res.empty();
-                }
+                createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
+                                               geodSrc, geodDst, context);
+                doFilterAndCheckPerfectOp = !res.empty();
             }
+        }
 
-            // NAD27 to NAD83 has tens of results already. No need to look
-            // for a pivot
-            if (!sameGeodeticDatum &&
-                ((res.empty() &&
-                  context.context->getAllowUseIntermediateCRS() ==
-                      CoordinateOperationContext::IntermediateCRSUse::
-                          IF_NO_DIRECT_TRANSFORMATION) ||
-                 context.context->getAllowUseIntermediateCRS() ==
-                     CoordinateOperationContext::IntermediateCRSUse::ALWAYS ||
-                 getenv("PROJ_FORCE_SEARCH_PIVOT"))) {
-                auto resWithIntermediate = findsOpsInRegistryWithIntermediate(
-                    sourceCRS, targetCRS, context.context);
-                res.insert(res.end(), resWithIntermediate.begin(),
-                           resWithIntermediate.end());
+        // NAD27 to NAD83 has tens of results already. No need to look
+        // for a pivot
+        if (!sameGeodeticDatum &&
+            ((res.empty() &&
+              context.context->getAllowUseIntermediateCRS() ==
+                  CoordinateOperationContext::IntermediateCRSUse::
+                      IF_NO_DIRECT_TRANSFORMATION) ||
+             context.context->getAllowUseIntermediateCRS() ==
+                 CoordinateOperationContext::IntermediateCRSUse::ALWAYS ||
+             getenv("PROJ_FORCE_SEARCH_PIVOT"))) {
+            auto resWithIntermediate = findsOpsInRegistryWithIntermediate(
+                sourceCRS, targetCRS, context.context);
+            res.insert(res.end(), resWithIntermediate.begin(),
+                       resWithIntermediate.end());
+            doFilterAndCheckPerfectOp = !res.empty();
+
+            // If transforming from/to WGS84 (Gxxxx), try through 'neutral'
+            // WGS84
+            if (res.empty() && geodSrc && geodDst &&
+                !context.inCreateOperationsThroughPreferredHub &&
+                !context.inCreateOperationsWithDatumPivotAntiRecursion) {
+                createOperationsThroughPreferredHub(res, sourceCRS, targetCRS,
+                                                    geodSrc, geodDst, context);
                 doFilterAndCheckPerfectOp = !res.empty();
+            }
+        }
+    }
+
+    if (doFilterAndCheckPerfectOp) {
+        // If we get at least a result with perfect accuracy, do not bother
+        // generating synthetic transforms.
+        if (hasPerfectAccuracyResult(res, context)) {
+            return true;
+        }
+    }
+    return false;
+}
 
-                // If transforming from/to WGS84 (Gxxxx), try through 'neutral'
-                // WGS84
-                if (res.empty() && geodSrc && geodDst &&
-                    !context.inCreateOperationsThroughPreferredHub &&
-                    !context.inCreateOperationsWithDatumPivotAntiRecursion) {
-                    createOperationsThroughPreferredHub(
-                        res, sourceCRS, targetCRS, geodSrc, geodDst, context);
-                    doFilterAndCheckPerfectOp = !res.empty();
+// ---------------------------------------------------------------------------
+
+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);
                 }
             }
         }
-
-        if (doFilterAndCheckPerfectOp) {
-            // If we get at least a result with perfect accuracy, do not bother
-            // generating synthetic transforms.
-            if (hasPerfectAccuracyResult(res, context)) {
-                return res;
+    } 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;
+}
 
-    // Special case if both CRS are geodetic
-    if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
+// ---------------------------------------------------------------------------
 
-        if (geodSrc->ellipsoid()->celestialBody() !=
-            geodDst->ellipsoid()->celestialBody()) {
-            throw util::UnsupportedOperationException(
-                "Source and target ellipsoid do not belong to the same "
-                "celestial body");
-        }
-
-        if (geogSrc && geogDst) {
-            return createOperationsGeogToGeog(res, sourceCRS, targetCRS,
-                                              geogSrc, geogDst);
-        }
-
-        const bool isSrcGeocentric = geodSrc->isGeocentric();
-        const bool isSrcGeographic = geogSrc != nullptr;
-        const bool isTargetGeocentric = geodDst->isGeocentric();
-        const bool isTargetGeographic = geogDst != nullptr;
-        if (((isSrcGeocentric && isTargetGeographic) ||
-             (isSrcGeographic && isTargetGeocentric)) &&
-            geodSrc->datum() != nullptr && geodDst->datum() != nullptr) {
-
-            // Same datum ?
-            if (geodSrc->datum()->_isEquivalentTo(
-                    geodDst->datum().get(),
-                    util::IComparable::Criterion::EQUIVALENT)) {
-                res.emplace_back(Conversion::createGeographicGeocentric(
-                    sourceCRS, targetCRS));
-            } else if (isSrcGeocentric) {
-                std::string interm_crs_name(geogDst->nameStr());
-                interm_crs_name += " (geocentric)";
-                auto interm_crs = util::nn_static_pointer_cast<crs::CRS>(
-                    crs::GeodeticCRS::create(
-                        addDomains(util::PropertyMap().set(
-                                       common::IdentifiedObject::NAME_KEY,
-                                       interm_crs_name),
-                                   geogDst),
-                        NN_NO_CHECK(geogDst->datum()),
-                        NN_CHECK_ASSERT(
-                            util::nn_dynamic_pointer_cast<cs::CartesianCS>(
-                                geodSrc->coordinateSystem()))));
-                auto opFirst =
-                    createBallparkGeocentricTranslation(sourceCRS, interm_crs);
-                auto opSecond = Conversion::createGeographicGeocentric(
-                    interm_crs, targetCRS);
-                res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-                    {opFirst, opSecond}, !allowEmptyIntersection));
-            } else {
-                return applyInverse(
-                    createOperations(targetCRS, sourceCRS, context));
+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;
+    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;
         }
 
-        if (isSrcGeocentric && isTargetGeocentric) {
-            res.emplace_back(
-                createBallparkGeocentricTranslation(sourceCRS, targetCRS));
-            return res;
+        ~AntiRecursionGuard() {
+            context.inCreateOperationsGeogToVertWithAlternativeGeog = false;
         }
+    };
+    AntiRecursionGuard guard(context);
 
-        // Transformation between two geodetic systems of unknown type
-        // This should normally not be triggered with "standard" CRS
-        res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS));
-        return res;
-    }
+    for (int i = 0; i < 2; i++) {
 
-    // If the source is a derived CRS, then chain the inverse of its
-    // deriving conversion, with transforms from its baseCRS to the
-    // targetCRS
-    if (derivedSrc) {
-        auto opFirst = derivedSrc->derivingConversion()->inverse();
-        // Small optimization if the targetCRS is the baseCRS of the source
-        // derivedCRS.
-        if (derivedSrc->baseCRS()->_isEquivalentTo(
-                targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) {
-            res.emplace_back(opFirst);
-            return res;
+        // 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());
+            }
         }
-        auto opsSecond =
-            createOperations(derivedSrc->baseCRS(), targetCRS, context);
-        for (const auto &opSecond : opsSecond) {
-            try {
-                res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-                    {opFirst, opSecond}, !allowEmptyIntersection));
-            } catch (const InvalidOperationEmptyIntersection &) {
+        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 res;
+        return false;
+    };
+
+    if (geog3DToVertTryThroughGeog2D(geogSrc, vertDst, targetCRS)) {
+        // do nothing
+    } else if (geog3DToVertTryThroughGeog2D(geogDst, vertSrc, sourceCRS)) {
+        res = applyInverse(res);
     }
 
-    // reverse of previous case
-    if (derivedDst) {
-        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    // 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 (boundSrc && geogDst) {
-        const auto &hubSrc = boundSrc->hubCRS();
-        auto hubSrcGeog =
-            dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
-        auto geogCRSOfBaseOfBoundSrc =
-            boundSrc->baseCRS()->extractGeographicCRS();
-        bool triedBoundCrsToGeogCRSSameAsHubCRS = false;
-        // Is it: boundCRS to a geogCRS that is the same as the hubCRS ?
-        if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
-            (hubSrcGeog->_isEquivalentTo(
-                 geogDst, util::IComparable::Criterion::EQUIVALENT) ||
-             hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) {
-            triedBoundCrsToGeogCRSSameAsHubCRS = true;
-            if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) {
-                // Optimization to avoid creating a useless concatenated
-                // operation
-                res.emplace_back(boundSrc->transformation());
-                return res;
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::GeodeticCRS *geodSrc,
+    const crs::GeodeticCRS *geodDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    if (geodSrc->ellipsoid()->celestialBody() !=
+        geodDst->ellipsoid()->celestialBody()) {
+        throw util::UnsupportedOperationException(
+            "Source and target ellipsoid do not belong to the same "
+            "celestial body");
+    }
+
+    auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(geodSrc);
+    auto geogDst = dynamic_cast<const crs::GeographicCRS *>(geodDst);
+
+    if (geogSrc && geogDst) {
+        createOperationsGeogToGeog(res, sourceCRS, targetCRS, geogSrc, geogDst);
+        return;
+    }
+
+    const bool isSrcGeocentric = geodSrc->isGeocentric();
+    const bool isSrcGeographic = geogSrc != nullptr;
+    const bool isTargetGeocentric = geodDst->isGeocentric();
+    const bool isTargetGeographic = geogDst != nullptr;
+    if (((isSrcGeocentric && isTargetGeographic) ||
+         (isSrcGeographic && isTargetGeocentric)) &&
+        geodSrc->datum() != nullptr && geodDst->datum() != nullptr) {
+
+        // Same datum ?
+        if (geodSrc->datum()->_isEquivalentTo(
+                geodDst->datum().get(),
+                util::IComparable::Criterion::EQUIVALENT)) {
+            res.emplace_back(
+                Conversion::createGeographicGeocentric(sourceCRS, targetCRS));
+        } else if (isSrcGeocentric) {
+            std::string interm_crs_name(geogDst->nameStr());
+            interm_crs_name += " (geocentric)";
+            auto interm_crs =
+                util::nn_static_pointer_cast<crs::CRS>(crs::GeodeticCRS::create(
+                    addDomains(util::PropertyMap().set(
+                                   common::IdentifiedObject::NAME_KEY,
+                                   interm_crs_name),
+                               geogDst),
+                    NN_NO_CHECK(geogDst->datum()),
+                    NN_CHECK_ASSERT(
+                        util::nn_dynamic_pointer_cast<cs::CartesianCS>(
+                            geodSrc->coordinateSystem()))));
+            auto opFirst =
+                createBallparkGeocentricTranslation(sourceCRS, interm_crs);
+            auto opSecond =
+                Conversion::createGeographicGeocentric(interm_crs, targetCRS);
+            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                {opFirst, opSecond}, !allowEmptyIntersection));
+        } else {
+            res = applyInverse(createOperations(targetCRS, sourceCRS, context));
+        }
+
+        return;
+    }
+
+    if (isSrcGeocentric && isTargetGeocentric) {
+        res.emplace_back(
+            createBallparkGeocentricTranslation(sourceCRS, targetCRS));
+        return;
+    }
+
+    // Transformation between two geodetic systems of unknown type
+    // This should normally not be triggered with "standard" CRS
+    res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS));
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsDerivedTo(
+    const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::DerivedCRS *derivedSrc,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    auto opFirst = derivedSrc->derivingConversion()->inverse();
+    // Small optimization if the targetCRS is the baseCRS of the source
+    // derivedCRS.
+    if (derivedSrc->baseCRS()->_isEquivalentTo(
+            targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) {
+        res.emplace_back(opFirst);
+        return;
+    }
+    auto opsSecond =
+        createOperations(derivedSrc->baseCRS(), targetCRS, context);
+    for (const auto &opSecond : opsSecond) {
+        try {
+            res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                {opFirst, opSecond}, !allowEmptyIntersection));
+        } catch (const InvalidOperationEmptyIntersection &) {
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::BoundCRS *boundSrc,
+    const crs::GeographicCRS *geogDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    const auto &hubSrc = boundSrc->hubCRS();
+    auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
+    auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS();
+    bool triedBoundCrsToGeogCRSSameAsHubCRS = false;
+    // Is it: boundCRS to a geogCRS that is the same as the hubCRS ?
+    if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
+        (hubSrcGeog->_isEquivalentTo(
+             geogDst, util::IComparable::Criterion::EQUIVALENT) ||
+         hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst)))) {
+        triedBoundCrsToGeogCRSSameAsHubCRS = true;
+        if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) {
+            // Optimization to avoid creating a useless concatenated
+            // operation
+            res.emplace_back(boundSrc->transformation());
+            return;
+        }
+        auto opsFirst = createOperations(
+            boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
+        if (!opsFirst.empty()) {
+            for (const auto &opFirst : opsFirst) {
+                try {
+                    res.emplace_back(
+                        ConcatenatedOperation::createComputeMetadata(
+                            {opFirst, boundSrc->transformation()},
+                            !allowEmptyIntersection));
+                } catch (const InvalidOperationEmptyIntersection &) {
+                }
             }
-            auto opsFirst =
-                createOperations(boundSrc->baseCRS(),
-                                 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
-            if (!opsFirst.empty()) {
-                for (const auto &opFirst : opsFirst) {
+            if (!res.empty()) {
+                return;
+            }
+        }
+        // If the datum are equivalent, this is also fine
+    } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+               geogDst->datum() &&
+               hubSrcGeog->datum()->_isEquivalentTo(
+                   geogDst->datum().get(),
+                   util::IComparable::Criterion::EQUIVALENT)) {
+        auto opsFirst = createOperations(
+            boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
+        auto opsLast = createOperations(hubSrc, targetCRS, context);
+        if (!opsFirst.empty() && !opsLast.empty()) {
+            for (const auto &opFirst : opsFirst) {
+                for (const auto &opLast : opsLast) {
                     try {
                         res.emplace_back(
                             ConcatenatedOperation::createComputeMetadata(
-                                {opFirst, boundSrc->transformation()},
+                                {opFirst, boundSrc->transformation(), opLast},
                                 !allowEmptyIntersection));
                     } catch (const InvalidOperationEmptyIntersection &) {
                     }
                 }
-                if (!res.empty()) {
-                    return res;
-                }
             }
-            // If the datum are equivalent, this is also fine
-        } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
-                   geogDst->datum() &&
-                   hubSrcGeog->datum()->_isEquivalentTo(
-                       geogDst->datum().get(),
-                       util::IComparable::Criterion::EQUIVALENT)) {
-            auto opsFirst =
-                createOperations(boundSrc->baseCRS(),
-                                 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
-            auto opsLast = createOperations(hubSrc, targetCRS, context);
-            if (!opsFirst.empty() && !opsLast.empty()) {
+            if (!res.empty()) {
+                return;
+            }
+        }
+        // Consider WGS 84 and NAD83 as equivalent in that context if the
+        // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27)
+        // Case of "+proj=latlong +ellps=clrk66
+        // +nadgrids=ntv1_can.dat,conus"
+        // to "+proj=latlong +datum=NAD83"
+    } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+               geogDst->datum() &&
+               geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
+                   datum::Ellipsoid::CLARKE_1866.get(),
+                   util::IComparable::Criterion::EQUIVALENT) &&
+               hubSrcGeog->datum()->_isEquivalentTo(
+                   datum::GeodeticReferenceFrame::EPSG_6326.get(),
+                   util::IComparable::Criterion::EQUIVALENT) &&
+               geogDst->datum()->_isEquivalentTo(
+                   datum::GeodeticReferenceFrame::EPSG_6269.get(),
+                   util::IComparable::Criterion::EQUIVALENT)) {
+        auto nnGeogCRSOfBaseOfBoundSrc = NN_NO_CHECK(geogCRSOfBaseOfBoundSrc);
+        if (boundSrc->baseCRS()->_isEquivalentTo(
+                nnGeogCRSOfBaseOfBoundSrc.get(),
+                util::IComparable::Criterion::EQUIVALENT)) {
+            auto transf = boundSrc->transformation()->shallowClone();
+            transf->setProperties(util::PropertyMap().set(
+                common::IdentifiedObject::NAME_KEY,
+                buildTransfName(boundSrc->baseCRS()->nameStr(),
+                                targetCRS->nameStr())));
+            transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr);
+            res.emplace_back(transf);
+            return;
+        } else {
+            auto opsFirst = createOperations(
+                boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context);
+            auto transf = boundSrc->transformation()->shallowClone();
+            transf->setProperties(util::PropertyMap().set(
+                common::IdentifiedObject::NAME_KEY,
+                buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(),
+                                targetCRS->nameStr())));
+            transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr);
+            if (!opsFirst.empty()) {
                 for (const auto &opFirst : opsFirst) {
-                    for (const auto &opLast : opsLast) {
-                        try {
-                            res.emplace_back(
-                                ConcatenatedOperation::createComputeMetadata(
-                                    {opFirst, boundSrc->transformation(),
-                                     opLast},
-                                    !allowEmptyIntersection));
-                        } catch (const InvalidOperationEmptyIntersection &) {
-                        }
+                    try {
+                        res.emplace_back(
+                            ConcatenatedOperation::createComputeMetadata(
+                                {opFirst, transf}, !allowEmptyIntersection));
+                    } catch (const InvalidOperationEmptyIntersection &) {
                     }
                 }
                 if (!res.empty()) {
-                    return res;
+                    return;
                 }
             }
-            // Consider WGS 84 and NAD83 as equivalent in that context if the
-            // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27)
-            // Case of "+proj=latlong +ellps=clrk66
-            // +nadgrids=ntv1_can.dat,conus"
-            // to "+proj=latlong +datum=NAD83"
-        } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
-                   geogDst->datum() &&
-                   geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
-                       datum::Ellipsoid::CLARKE_1866.get(),
-                       util::IComparable::Criterion::EQUIVALENT) &&
-                   hubSrcGeog->datum()->_isEquivalentTo(
-                       datum::GeodeticReferenceFrame::EPSG_6326.get(),
-                       util::IComparable::Criterion::EQUIVALENT) &&
-                   geogDst->datum()->_isEquivalentTo(
-                       datum::GeodeticReferenceFrame::EPSG_6269.get(),
-                       util::IComparable::Criterion::EQUIVALENT)) {
-            auto nnGeogCRSOfBaseOfBoundSrc =
-                NN_NO_CHECK(geogCRSOfBaseOfBoundSrc);
-            if (boundSrc->baseCRS()->_isEquivalentTo(
-                    nnGeogCRSOfBaseOfBoundSrc.get(),
-                    util::IComparable::Criterion::EQUIVALENT)) {
-                auto transf = boundSrc->transformation()->shallowClone();
-                transf->setProperties(util::PropertyMap().set(
-                    common::IdentifiedObject::NAME_KEY,
-                    buildTransfName(boundSrc->baseCRS()->nameStr(),
-                                    targetCRS->nameStr())));
-                transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr);
-                res.emplace_back(transf);
-                return res;
-            } else {
-                auto opsFirst = createOperations(
-                    boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context);
-                auto transf = boundSrc->transformation()->shallowClone();
-                transf->setProperties(util::PropertyMap().set(
-                    common::IdentifiedObject::NAME_KEY,
-                    buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(),
-                                    targetCRS->nameStr())));
-                transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr);
-                if (!opsFirst.empty()) {
-                    for (const auto &opFirst : opsFirst) {
+        }
+    }
+
+    if (hubSrcGeog &&
+        hubSrcGeog->_isEquivalentTo(geogDst,
+                                    util::IComparable::Criterion::EQUIVALENT) &&
+        dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get())) {
+        res.emplace_back(boundSrc->transformation());
+        return;
+    }
+
+    if (!triedBoundCrsToGeogCRSSameAsHubCRS && hubSrcGeog &&
+        geogCRSOfBaseOfBoundSrc) {
+        // This one should go to the above 'Is it: boundCRS to a geogCRS
+        // that is the same as the hubCRS ?' case
+        auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+        auto opsLast = createOperations(hubSrc, targetCRS, context);
+        if (!opsFirst.empty() && !opsLast.empty()) {
+            for (const auto &opFirst : opsFirst) {
+                for (const auto &opLast : opsLast) {
+                    // Exclude artificial transformations from the hub
+                    // to the target CRS
+                    if (!opLast->hasBallparkTransformation()) {
                         try {
                             res.emplace_back(
                                 ConcatenatedOperation::createComputeMetadata(
-                                    {opFirst, transf},
+                                    {opFirst, opLast},
                                     !allowEmptyIntersection));
                         } catch (const InvalidOperationEmptyIntersection &) {
                         }
                     }
-                    if (!res.empty()) {
-                        return res;
-                    }
                 }
             }
+            if (!res.empty()) {
+                return;
+            }
         }
+    }
 
-        if (hubSrcGeog &&
-            hubSrcGeog->_isEquivalentTo(
-                geogDst, util::IComparable::Criterion::EQUIVALENT) &&
-            dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get())) {
-            res.emplace_back(boundSrc->transformation());
-            return res;
-        }
-
-        if (!triedBoundCrsToGeogCRSSameAsHubCRS && hubSrcGeog &&
-            geogCRSOfBaseOfBoundSrc) {
-            // This one should go to the above 'Is it: boundCRS to a geogCRS
-            // that is the same as the hubCRS ?' case
-            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
-            auto opsLast = createOperations(hubSrc, targetCRS, context);
-            if (!opsFirst.empty() && !opsLast.empty()) {
+    auto vertCRSOfBaseOfBoundSrc =
+        dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
+    if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) {
+        auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+        if (context.skipHorizontalTransformation) {
+            if (!opsFirst.empty())
+                res = opsFirst;
+            return;
+        } else {
+            auto opsSecond = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsSecond.empty()) {
                 for (const auto &opFirst : opsFirst) {
-                    for (const auto &opLast : opsLast) {
+                    for (const auto &opLast : opsSecond) {
                         // Exclude artificial transformations from the hub
                         // to the target CRS
                         if (!opLast->hasBallparkTransformation()) {
@@ -12965,651 +13315,728 @@ CoordinateOperationFactory::Private::createOperations(
                     }
                 }
                 if (!res.empty()) {
-                    return res;
+                    return;
                 }
             }
         }
+    }
 
-        auto vertCRSOfBaseOfBoundSrc =
-            dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
-        if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) {
-            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
-            if (context.skipHorizontalTransformation) {
-                if (!opsFirst.empty())
-                    return opsFirst;
-            } else {
-                auto opsSecond = createOperations(hubSrc, targetCRS, context);
-                if (!opsFirst.empty() && !opsSecond.empty()) {
-                    for (const auto &opFirst : opsFirst) {
-                        for (const auto &opLast : opsSecond) {
-                            // Exclude artificial transformations from the hub
-                            // to the target CRS
-                            if (!opLast->hasBallparkTransformation()) {
-                                try {
-                                    res.emplace_back(
-                                        ConcatenatedOperation::
-                                            createComputeMetadata(
-                                                {opFirst, opLast},
-                                                !allowEmptyIntersection));
-                                } catch (
-                                    const InvalidOperationEmptyIntersection &) {
-                                }
-                            }
-                        }
-                    }
-                    if (!res.empty()) {
-                        return res;
-                    }
-                }
-            }
-        }
+    res = createOperations(boundSrc->baseCRS(), targetCRS, context);
+}
 
-        return createOperations(boundSrc->baseCRS(), targetCRS, context);
-    }
+// ---------------------------------------------------------------------------
 
-    // reverse of previous case
-    if (geogSrc && boundDst) {
-        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+void CoordinateOperationFactory::Private::createOperationsBoundToVert(
+    const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::BoundCRS *boundSrc,
+    const crs::VerticalCRS *vertDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    auto baseSrcVert =
+        dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
+    const auto &hubSrc = boundSrc->hubCRS();
+    auto hubSrcVert = dynamic_cast<const crs::VerticalCRS *>(hubSrc.get());
+    if (baseSrcVert && hubSrcVert &&
+        vertDst->_isEquivalentTo(hubSrcVert,
+                                 util::IComparable::Criterion::EQUIVALENT)) {
+        res.emplace_back(boundSrc->transformation());
+        return;
     }
 
-    // vertCRS (as boundCRS with transformation to target vertCRS) to
-    // vertCRS
-    if (boundSrc && vertDst) {
-        auto baseSrcVert =
-            dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
-        const auto &hubSrc = boundSrc->hubCRS();
-        auto hubSrcVert = dynamic_cast<const crs::VerticalCRS *>(hubSrc.get());
-        if (baseSrcVert && hubSrcVert &&
-            vertDst->_isEquivalentTo(
-                hubSrcVert, util::IComparable::Criterion::EQUIVALENT)) {
-            res.emplace_back(boundSrc->transformation());
-            return res;
-        }
+    res = createOperations(boundSrc->baseCRS(), targetCRS, context);
+}
 
-        return createOperations(boundSrc->baseCRS(), targetCRS, context);
-    }
+// ---------------------------------------------------------------------------
 
-    // reverse of previous case
-    if (boundDst && vertSrc) {
-        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+void CoordinateOperationFactory::Private::createOperationsVertToVert(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context & /*context*/, const crs::VerticalCRS *vertSrc,
+    const crs::VerticalCRS *vertDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    const auto &srcDatum = vertSrc->datum();
+    const auto &dstDatum = vertDst->datum();
+    const bool equivalentVDatum =
+        (srcDatum && dstDatum &&
+         srcDatum->_isEquivalentTo(dstDatum.get(),
+                                   util::IComparable::Criterion::EQUIVALENT));
+
+    const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0];
+    const double convSrc = srcAxis->unit().conversionToSI();
+    const auto &dstAxis = vertDst->coordinateSystem()->axisList()[0];
+    const double convDst = dstAxis->unit().conversionToSI();
+    const bool srcIsUp = srcAxis->direction() == cs::AxisDirection::UP;
+    const bool srcIsDown = srcAxis->direction() == cs::AxisDirection::DOWN;
+    const bool dstIsUp = dstAxis->direction() == cs::AxisDirection::UP;
+    const bool dstIsDown = dstAxis->direction() == cs::AxisDirection::DOWN;
+    const bool heightDepthReversal =
+        ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp));
+
+    const double factor = convSrc / convDst;
+    auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
+    if (!equivalentVDatum) {
+        name += BALLPARK_VERTICAL_TRANSFORMATION;
+        auto conv = Transformation::createChangeVerticalUnit(
+            util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name),
+            sourceCRS, targetCRS,
+            // In case of a height depth reversal, we should probably have
+            // 2 steps instead of putting a negative factor...
+            common::Scale(heightDepthReversal ? -factor : factor), {});
+        conv->setHasBallparkTransformation(true);
+        res.push_back(conv);
+    } else if (convSrc != convDst) {
+        auto conv = Conversion::createChangeVerticalUnit(
+            util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name),
+            // In case of a height depth reversal, we should probably have
+            // 2 steps instead of putting a negative factor...
+            common::Scale(heightDepthReversal ? -factor : factor));
+        conv->setCRSs(sourceCRS, targetCRS, nullptr);
+        res.push_back(conv);
+    } else if (heightDepthReversal) {
+        auto conv = Conversion::createHeightDepthReversal(
+            util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name));
+        conv->setCRSs(sourceCRS, targetCRS, nullptr);
+        res.push_back(conv);
     }
+}
 
-    if (vertSrc && vertDst) {
-        const auto &srcDatum = vertSrc->datum();
-        const auto &dstDatum = vertDst->datum();
-        const bool equivalentVDatum =
-            (srcDatum && dstDatum &&
-             srcDatum->_isEquivalentTo(
-                 dstDatum.get(), util::IComparable::Criterion::EQUIVALENT));
-
-        const double convSrc =
-            vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
-        const double convDst =
-            vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI();
-
-        const double factor = convSrc / convDst;
-        auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
-        if (!equivalentVDatum) {
-            name += BALLPARK_VERTICAL_TRANSFORMATION;
-            auto conv = Transformation::createChangeVerticalUnit(
-                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                        name),
-                sourceCRS, targetCRS, common::Scale(factor), {});
-            conv->setHasBallparkTransformation(true);
-            res.push_back(conv);
-        } else if (convSrc != convDst) {
-            auto conv = Conversion::createChangeVerticalUnit(
-                util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
-                                        name),
-                common::Scale(factor));
-            conv->setCRSs(sourceCRS, targetCRS, nullptr);
-            res.push_back(conv);
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsVertToGeog(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::VerticalCRS *vertSrc,
+    const crs::GeographicCRS *geogDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    if (vertSrc->identifiers().empty()) {
+        const auto &vertSrcName = vertSrc->nameStr();
+        const auto &authFactory = context.context->getAuthorityFactory();
+        if (authFactory != nullptr && vertSrcName != "unnamed" &&
+            vertSrcName != "unknown") {
+            auto matches = authFactory->createObjectsFromName(
+                vertSrcName, {io::AuthorityFactory::ObjectType::VERTICAL_CRS},
+                false, 2);
+            if (matches.size() == 1) {
+                const auto &match = matches.front();
+                if (vertSrc->_isEquivalentTo(
+                        match.get(),
+                        util::IComparable::Criterion::EQUIVALENT) &&
+                    !match->identifiers().empty()) {
+                    res = createOperations(
+                        NN_NO_CHECK(
+                            util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
+                                match)),
+                        targetCRS, context);
+                    return;
+                }
+            }
         }
-        return res;
     }
 
-    // A bit odd case as we are comparing apples to oranges, but in case
-    // the vertical unit differ, do something useful.
-    if (vertSrc && geogDst) {
+    const double convSrc =
+        vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+    double convDst = 1.0;
+    const auto &geogAxis = geogDst->coordinateSystem()->axisList();
+    if (geogAxis.size() == 3) {
+        convDst = geogAxis[2]->unit().conversionToSI();
+    }
 
-        if (vertSrc->identifiers().empty()) {
-            const auto &vertSrcName = vertSrc->nameStr();
-            if (authFactory != nullptr && vertSrcName != "unnamed" &&
-                vertSrcName != "unknown") {
-                auto matches = authFactory->createObjectsFromName(
-                    vertSrcName,
-                    {io::AuthorityFactory::ObjectType::VERTICAL_CRS}, false, 2);
-                if (matches.size() == 1) {
-                    const auto &match = matches.front();
-                    if (vertSrc->_isEquivalentTo(
-                            match.get(),
-                            util::IComparable::Criterion::EQUIVALENT) &&
-                        !match->identifiers().empty()) {
-                        return createOperations(
-                            NN_NO_CHECK(
-                                util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
-                                    match)),
-                            targetCRS, context);
+    const double factor = convSrc / convDst;
+    auto conv = Transformation::createChangeVerticalUnit(
+        util::PropertyMap().set(
+            common::IdentifiedObject::NAME_KEY,
+            buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
+                BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
+        sourceCRS, targetCRS, common::Scale(factor), {});
+    conv->setHasBallparkTransformation(true);
+    res.push_back(conv);
+}
+
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsBoundToBound(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::BoundCRS *boundSrc,
+    const crs::BoundCRS *boundDst, std::vector<CoordinateOperationNNPtr> &res) {
+
+    const auto &hubSrc = boundSrc->hubCRS();
+    auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
+    const auto &hubDst = boundDst->hubCRS();
+    auto hubDstGeog = dynamic_cast<const crs::GeographicCRS *>(hubDst.get());
+    auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS();
+    auto geogCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractGeographicCRS();
+    if (hubSrcGeog && hubDstGeog &&
+        hubSrcGeog->_isEquivalentTo(hubDstGeog,
+                                    util::IComparable::Criterion::EQUIVALENT) &&
+        geogCRSOfBaseOfBoundSrc && geogCRSOfBaseOfBoundDst) {
+        const bool firstIsNoOp = geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
+            boundSrc->baseCRS().get(),
+            util::IComparable::Criterion::EQUIVALENT);
+        const bool lastIsNoOp = geogCRSOfBaseOfBoundDst->_isEquivalentTo(
+            boundDst->baseCRS().get(),
+            util::IComparable::Criterion::EQUIVALENT);
+        auto opsFirst = createOperations(
+            boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
+        auto opsLast = createOperations(NN_NO_CHECK(geogCRSOfBaseOfBoundDst),
+                                        boundDst->baseCRS(), context);
+        if (!opsFirst.empty() && !opsLast.empty()) {
+            const auto &opSecond = boundSrc->transformation();
+            auto opThird = boundDst->transformation()->inverse();
+            for (const auto &opFirst : opsFirst) {
+                for (const auto &opLast : opsLast) {
+                    try {
+                        std::vector<CoordinateOperationNNPtr> ops;
+                        if (!firstIsNoOp) {
+                            ops.push_back(opFirst);
+                        }
+                        ops.push_back(opSecond);
+                        ops.push_back(opThird);
+                        if (!lastIsNoOp) {
+                            ops.push_back(opLast);
+                        }
+                        res.emplace_back(
+                            ConcatenatedOperation::createComputeMetadata(
+                                ops, !allowEmptyIntersection));
+                    } catch (const InvalidOperationEmptyIntersection &) {
                     }
                 }
             }
+            if (!res.empty()) {
+                return;
+            }
         }
+    }
 
-        const double convSrc =
-            vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
-        double convDst = 1.0;
-        const auto &geogAxis = geogDst->coordinateSystem()->axisList();
-        if (geogAxis.size() == 3) {
-            convDst = geogAxis[2]->unit().conversionToSI();
+    auto vertCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractVerticalCRS();
+    auto vertCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractVerticalCRS();
+    if (hubSrcGeog && hubDstGeog &&
+        hubSrcGeog->_isEquivalentTo(hubDstGeog,
+                                    util::IComparable::Criterion::EQUIVALENT) &&
+        vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) {
+        auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+        auto opsLast = createOperations(hubSrc, targetCRS, context);
+        if (!opsFirst.empty() && !opsLast.empty()) {
+            for (const auto &opFirst : opsFirst) {
+                for (const auto &opLast : opsLast) {
+                    try {
+                        res.emplace_back(
+                            ConcatenatedOperation::createComputeMetadata(
+                                {opFirst, opLast}, !allowEmptyIntersection));
+                    } catch (const InvalidOperationEmptyIntersection &) {
+                    }
+                }
+            }
+            if (!res.empty()) {
+                return;
+            }
         }
-
-        const double factor = convSrc / convDst;
-        auto conv = Transformation::createChangeVerticalUnit(
-            util::PropertyMap().set(
-                common::IdentifiedObject::NAME_KEY,
-                buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
-                    BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
-            sourceCRS, targetCRS, common::Scale(factor), {});
-        conv->setHasBallparkTransformation(true);
-        res.push_back(conv);
-        return res;
     }
 
-    // reverse of previous case
-    if (vertDst && geogSrc) {
-        return applyInverse(createOperations(targetCRS, sourceCRS, context));
-    }
+    res = createOperations(boundSrc->baseCRS(), boundDst->baseCRS(), context);
+}
 
-    // boundCRS to boundCRS using the same geographic hubCRS
-    if (boundSrc && boundDst) {
-        const auto &hubSrc = boundSrc->hubCRS();
-        auto hubSrcGeog =
-            dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
-        const auto &hubDst = boundDst->hubCRS();
-        auto hubDstGeog =
-            dynamic_cast<const crs::GeographicCRS *>(hubDst.get());
-        auto geogCRSOfBaseOfBoundSrc =
-            boundSrc->baseCRS()->extractGeographicCRS();
-        auto geogCRSOfBaseOfBoundDst =
-            boundDst->baseCRS()->extractGeographicCRS();
-        if (hubSrcGeog && hubDstGeog &&
-            hubSrcGeog->_isEquivalentTo(
-                hubDstGeog, util::IComparable::Criterion::EQUIVALENT) &&
-            geogCRSOfBaseOfBoundSrc && geogCRSOfBaseOfBoundDst) {
-            const bool firstIsNoOp = geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
-                boundSrc->baseCRS().get(),
-                util::IComparable::Criterion::EQUIVALENT);
-            const bool lastIsNoOp = geogCRSOfBaseOfBoundDst->_isEquivalentTo(
-                boundDst->baseCRS().get(),
-                util::IComparable::Criterion::EQUIVALENT);
-            auto opsFirst =
-                createOperations(boundSrc->baseCRS(),
-                                 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
-            auto opsLast =
-                createOperations(NN_NO_CHECK(geogCRSOfBaseOfBoundDst),
-                                 boundDst->baseCRS(), context);
-            if (!opsFirst.empty() && !opsLast.empty()) {
-                const auto &opSecond = boundSrc->transformation();
-                auto opThird = boundDst->transformation()->inverse();
-                for (const auto &opFirst : opsFirst) {
-                    for (const auto &opLast : opsLast) {
-                        try {
-                            std::vector<CoordinateOperationNNPtr> ops;
-                            if (!firstIsNoOp) {
-                                ops.push_back(opFirst);
-                            }
-                            ops.push_back(opSecond);
-                            ops.push_back(opThird);
-                            if (!lastIsNoOp) {
-                                ops.push_back(opLast);
-                            }
-                            res.emplace_back(
-                                ConcatenatedOperation::createComputeMetadata(
-                                    ops, !allowEmptyIntersection));
-                        } catch (const InvalidOperationEmptyIntersection &) {
-                        }
-                    }
+// ---------------------------------------------------------------------------
+
+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 (!res.empty()) {
-                    return res;
+
+                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;
+}
 
-        auto vertCRSOfBaseOfBoundSrc =
-            boundSrc->baseCRS()->extractVerticalCRS();
-        auto vertCRSOfBaseOfBoundDst =
-            boundDst->baseCRS()->extractVerticalCRS();
-        if (hubSrcGeog && hubDstGeog &&
-            hubSrcGeog->_isEquivalentTo(
-                hubDstGeog, util::IComparable::Criterion::EQUIVALENT) &&
-            vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) {
-            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
-            auto opsLast = createOperations(hubSrc, targetCRS, context);
-            if (!opsFirst.empty() && !opsLast.empty()) {
-                for (const auto &opFirst : opsFirst) {
-                    for (const auto &opLast : opsLast) {
-                        try {
-                            res.emplace_back(
-                                ConcatenatedOperation::createComputeMetadata(
-                                    {opFirst, opLast},
-                                    !allowEmptyIntersection));
-                        } catch (const InvalidOperationEmptyIntersection &) {
-                        }
+// ---------------------------------------------------------------------------
+
+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 and/or height depth reversal +
+            // transformation ?
+            for (const auto &step : steps) {
+                const auto transf =
+                    dynamic_cast<const Transformation *>(step.get());
+                if (transf) {
+                    // Only support a single Transformation in the steps
+                    if (transformationVerticalTransform != nullptr) {
+                        transformationVerticalTransform = nullptr;
+                        break;
                     }
-                }
-                if (!res.empty()) {
-                    return res;
+                    transformationVerticalTransform = transf;
                 }
             }
         }
+    }
+    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());
+        }
+    }
 
-        return createOperations(boundSrc->baseCRS(), boundDst->baseCRS(),
-                                context);
+    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();
+        }
     }
 
-    auto compoundSrc = dynamic_cast<crs::CompoundCRS *>(sourceCRS.get());
-    if (compoundSrc && geogDst) {
-        const auto &componentsSrc = compoundSrc->componentReferenceSystems();
-        if (!componentsSrc.empty()) {
-            std::vector<CoordinateOperationNNPtr> horizTransforms;
-            auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
-            if (srcGeogCRS) {
-                horizTransforms =
-                    createOperations(componentsSrc[0], targetCRS, context);
-            }
-            std::vector<CoordinateOperationNNPtr> verticalTransforms;
-            if (componentsSrc.size() >= 2 &&
-                componentsSrc[1]->extractVerticalCRS()) {
+    return interpolationGeogCRS;
+}
 
-                struct SetSkipHorizontalTransform {
-                    Context &context;
+// ---------------------------------------------------------------------------
 
-                    explicit SetSkipHorizontalTransform(Context &contextIn)
-                        : context(contextIn) {
-                        assert(!context.skipHorizontalTransformation);
-                        context.skipHorizontalTransformation = true;
-                    }
+void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::CompoundCRS *compoundSrc,
+    const crs::GeographicCRS *geogDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
 
-                    ~SetSkipHorizontalTransform() {
-                        context.skipHorizontalTransformation = false;
-                    }
-                };
-                SetSkipHorizontalTransform setSkipHorizontalTransform(context);
-
-                verticalTransforms =
-                    createOperations(componentsSrc[1], targetCRS, context);
-                bool foundRegisteredTransformWithAllGridsAvailable = false;
-                for (const auto &op : verticalTransforms) {
-                    if (!op->identifiers().empty() && authFactory) {
-                        bool missingGrid = false;
-                        const auto gridsNeeded =
-                            op->gridsNeeded(authFactory->databaseContext());
+    const auto &authFactory = context.context->getAuthorityFactory();
+    const auto &componentsSrc = compoundSrc->componentReferenceSystems();
+    if (!componentsSrc.empty()) {
+        std::vector<CoordinateOperationNNPtr> horizTransforms;
+        auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
+        if (srcGeogCRS) {
+            horizTransforms =
+                createOperations(componentsSrc[0], targetCRS, context);
+        }
+        std::vector<CoordinateOperationNNPtr> verticalTransforms;
+
+        const auto dbContext =
+            authFactory ? authFactory->databaseContext().as_nullable()
+                        : nullptr;
+        if (componentsSrc.size() >= 2 &&
+            componentsSrc[1]->extractVerticalCRS()) {
+
+            struct SetSkipHorizontalTransform {
+                Context &context;
+
+                explicit SetSkipHorizontalTransform(Context &contextIn)
+                    : context(contextIn) {
+                    assert(!context.skipHorizontalTransformation);
+                    context.skipHorizontalTransformation = true;
+                }
+
+                ~SetSkipHorizontalTransform() {
+                    context.skipHorizontalTransformation = false;
+                }
+            };
+            SetSkipHorizontalTransform setSkipHorizontalTransform(context);
+
+            verticalTransforms = createOperations(
+                componentsSrc[1],
+                targetCRS->promoteTo3D(std::string(), dbContext), context);
+            bool foundRegisteredTransformWithAllGridsAvailable = false;
+            const bool ignoreMissingGrids =
+                context.context->getGridAvailabilityUse() ==
+                CoordinateOperationContext::GridAvailabilityUse::
+                    IGNORE_GRID_AVAILABILITY;
+            for (const auto &op : verticalTransforms) {
+                if (hasIdentifiers(op) && dbContext) {
+                    bool missingGrid = false;
+                    if (!ignoreMissingGrids) {
+                        const auto gridsNeeded = op->gridsNeeded(dbContext);
                         for (const auto &gridDesc : gridsNeeded) {
                             if (!gridDesc.available) {
                                 missingGrid = true;
                                 break;
                             }
                         }
-                        if (!missingGrid) {
-                            foundRegisteredTransformWithAllGridsAvailable =
-                                true;
-                            break;
-                        }
+                    }
+                    if (!missingGrid) {
+                        foundRegisteredTransformWithAllGridsAvailable = true;
+                        break;
                     }
                 }
-                if (!foundRegisteredTransformWithAllGridsAvailable &&
-                    srcGeogCRS &&
-                    !srcGeogCRS->_isEquivalentTo(
-                        geogDst, util::IComparable::Criterion::EQUIVALENT) &&
-                    !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) {
-                    auto verticalTransformsTmp = createOperations(
-                        componentsSrc[1], NN_NO_CHECK(srcGeogCRS), context);
-                    bool foundRegisteredTransform = false;
-                    foundRegisteredTransformWithAllGridsAvailable = false;
-                    for (const auto &op : verticalTransformsTmp) {
-                        if (!op->identifiers().empty() && authFactory) {
-                            bool missingGrid = false;
-                            const auto gridsNeeded =
-                                op->gridsNeeded(authFactory->databaseContext());
+            }
+            if (!foundRegisteredTransformWithAllGridsAvailable && srcGeogCRS &&
+                !srcGeogCRS->_isEquivalentTo(
+                    geogDst, util::IComparable::Criterion::EQUIVALENT) &&
+                !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst))) {
+                auto verticalTransformsTmp = createOperations(
+                    componentsSrc[1],
+                    NN_NO_CHECK(srcGeogCRS)
+                        ->promoteTo3D(std::string(), dbContext),
+                    context);
+                bool foundRegisteredTransform = false;
+                foundRegisteredTransformWithAllGridsAvailable = false;
+                for (const auto &op : verticalTransformsTmp) {
+                    if (hasIdentifiers(op) && dbContext) {
+                        bool missingGrid = false;
+                        if (!ignoreMissingGrids) {
+                            const auto gridsNeeded = op->gridsNeeded(dbContext);
                             for (const auto &gridDesc : gridsNeeded) {
                                 if (!gridDesc.available) {
                                     missingGrid = true;
                                     break;
                                 }
                             }
-                            foundRegisteredTransform = true;
-                            if (!missingGrid) {
-                                foundRegisteredTransformWithAllGridsAvailable =
-                                    true;
-                                break;
-                            }
                         }
-                    }
-                    if (foundRegisteredTransformWithAllGridsAvailable) {
-                        verticalTransforms = verticalTransformsTmp;
-                    } else if (foundRegisteredTransform) {
-                        verticalTransforms.insert(verticalTransforms.end(),
-                                                  verticalTransformsTmp.begin(),
-                                                  verticalTransformsTmp.end());
+                        foundRegisteredTransform = true;
+                        if (!missingGrid) {
+                            foundRegisteredTransformWithAllGridsAvailable =
+                                true;
+                            break;
+                        }
                     }
                 }
+                if (foundRegisteredTransformWithAllGridsAvailable) {
+                    verticalTransforms = verticalTransformsTmp;
+                } else if (foundRegisteredTransform) {
+                    verticalTransforms.insert(verticalTransforms.end(),
+                                              verticalTransformsTmp.begin(),
+                                              verticalTransformsTmp.end());
+                }
             }
+        }
 
-            if (horizTransforms.empty() || verticalTransforms.empty()) {
-                return horizTransforms;
-            }
-
-            typedef std::pair<std::vector<CoordinateOperationNNPtr>,
-                              std::vector<CoordinateOperationNNPtr>>
-                PairOfTransforms;
-            std::map<std::string, PairOfTransforms>
-                cacheHorizToInterpAndInterpToTarget;
+        if (horizTransforms.empty() || verticalTransforms.empty()) {
+            res = horizTransforms;
+            return;
+        }
 
-            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());
-                    }
+        typedef std::pair<std::vector<CoordinateOperationNNPtr>,
+                          std::vector<CoordinateOperationNNPtr>>
+            PairOfTransforms;
+        std::map<std::string, PairOfTransforms>
+            cacheHorizToInterpAndInterpToTarget;
+
+        for (const auto &verticalTransform : verticalTransforms) {
+            crs::GeographicCRSPtr interpolationGeogCRS =
+                getInterpolationGeogCRS(verticalTransform, dbContext);
+            if (interpolationGeogCRS) {
+                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 (interpolationGeogCRS) {
-                    if (interpolationGeogCRS->coordinateSystem()
-                            ->axisList()
-                            .size() == 3) {
-                        io::DatabaseContextPtr dbContext;
-                        if (authFactory) {
-                            dbContext = authFactory->databaseContext();
-                        }
-                        // 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);
-                            interpToTargetOps = createOperations(
-                                NN_NO_CHECK(interpolationGeogCRS), targetCRS,
-                                context);
-                            cacheHorizToInterpAndInterpToTarget[key] =
-                                PairOfTransforms(srcToInterpOps,
-                                                 interpToTargetOps);
-                        } else {
-                            srcToInterpOps = iter->second.first;
-                            interpToTargetOps = iter->second.second;
-                        }
-                    } else {
+                if (!key.empty()) {
+                    auto iter = cacheHorizToInterpAndInterpToTarget.find(key);
+                    if (iter == cacheHorizToInterpAndInterpToTarget.end()) {
                         srcToInterpOps = createOperations(
                             componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS),
                             context);
-                        interpToTargetOps =
-                            createOperations(NN_NO_CHECK(interpolationGeogCRS),
-                                             targetCRS, context);
+                        interpToTargetOps = createOperations(
+                            NN_NO_CHECK(interpolationGeogCRS),
+                            targetCRS->demoteTo2D(std::string(), dbContext),
+                            context);
+                        cacheHorizToInterpAndInterpToTarget[key] =
+                            PairOfTransforms(srcToInterpOps, interpToTargetOps);
+                    } else {
+                        srcToInterpOps = iter->second.first;
+                        interpToTargetOps = iter->second.second;
                     }
+                } else {
+                    srcToInterpOps = createOperations(
+                        componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS),
+                        context);
+                    interpToTargetOps = createOperations(
+                        NN_NO_CHECK(interpolationGeogCRS),
+                        targetCRS->demoteTo2D(std::string(), dbContext),
+                        context);
+                }
 
-                    for (const auto &srcToInterp : srcToInterpOps) {
-                        for (const auto &interpToTarget : interpToTargetOps) {
-
-                            if (useDifferentTransformationsForSameSourceTarget(
-                                    srcToInterp, interpToTarget)) {
-                                continue;
-                            }
+                for (const auto &srcToInterp : srcToInterpOps) {
+                    for (const auto &interpToTarget : interpToTargetOps) {
 
-                            try {
-                                auto op = createHorizVerticalHorizPROJBased(
-                                    sourceCRS, targetCRS, srcToInterp,
-                                    verticalTransform, interpToTarget,
-                                    interpolationGeogCRS, true);
-                                res.emplace_back(op);
-                            } catch (const std::exception &) {
-                            }
+                        if (useDifferentTransformationsForSameSourceTarget(
+                                srcToInterp, interpToTarget)) {
+                            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) {
+
                         try {
-                            auto op = createHorizVerticalPROJBased(
-                                sourceCRS, targetCRS, horizTransform,
-                                verticalTransform);
+                            auto op = createHorizVerticalHorizPROJBased(
+                                sourceCRS, targetCRS, srcToInterp,
+                                verticalTransform, interpToTarget,
+                                interpolationGeogCRS, true);
                             res.emplace_back(op);
                         } catch (const std::exception &) {
                         }
                     }
                 }
-            }
-
-            return res;
-        }
-    } else if (compoundSrc && geodDst) {
-        auto datum = geodDst->datum();
-        if (datum) {
-            auto cs =
-                cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight(
-                    common::UnitOfMeasure::DEGREE,
-                    common::UnitOfMeasure::METRE);
-            auto intermGeog3DCRS = util::nn_static_pointer_cast<crs::CRS>(
-                crs::GeographicCRS::create(
-                    util::PropertyMap()
-                        .set(common::IdentifiedObject::NAME_KEY,
-                             geodDst->nameStr())
-                        .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
-                             metadata::Extent::WORLD),
-                    NN_NO_CHECK(datum), cs));
-            auto sourceToGeog3DOps =
-                createOperations(sourceCRS, intermGeog3DCRS, context);
-            auto geog3DToTargetOps =
-                createOperations(intermGeog3DCRS, targetCRS, context);
-            if (!geog3DToTargetOps.empty()) {
-                for (const auto &op : sourceToGeog3DOps) {
-                    res.emplace_back(
-                        ConcatenatedOperation::createComputeMetadata(
-                            {op, geog3DToTargetOps.front()},
-                            !allowEmptyIntersection));
+            } 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) {
+                    try {
+                        auto op = createHorizVerticalPROJBased(
+                            sourceCRS, targetCRS, horizTransform,
+                            verticalTransform);
+                        res.emplace_back(op);
+                    } catch (const std::exception &) {
+                    }
                 }
-                return res;
             }
         }
     }
+}
 
-    // reverse of previous case
-    auto compoundDst = dynamic_cast<const crs::CompoundCRS *>(targetCRS.get());
-    if (geodSrc && compoundDst) {
-        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsCompoundToGeod(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::CompoundCRS * /*compoundSrc*/,
+    const crs::GeodeticCRS *geodDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+    auto datum = geodDst->datum();
+    if (datum) {
+        auto cs = cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight(
+            common::UnitOfMeasure::DEGREE, common::UnitOfMeasure::METRE);
+        auto intermGeog3DCRS =
+            util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
+                util::PropertyMap()
+                    .set(common::IdentifiedObject::NAME_KEY, geodDst->nameStr())
+                    .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                         metadata::Extent::WORLD),
+                NN_NO_CHECK(datum), cs));
+        auto sourceToGeog3DOps =
+            createOperations(sourceCRS, intermGeog3DCRS, context);
+        auto geog3DToTargetOps =
+            createOperations(intermGeog3DCRS, targetCRS, context);
+        if (!geog3DToTargetOps.empty()) {
+            for (const auto &op : sourceToGeog3DOps) {
+                res.emplace_back(ConcatenatedOperation::createComputeMetadata(
+                    {op, geog3DToTargetOps.front()}, !allowEmptyIntersection));
+            }
+        }
     }
+}
 
-    if (compoundSrc && compoundDst) {
-        const auto &componentsSrc = compoundSrc->componentReferenceSystems();
-        const auto &componentsDst = compoundDst->componentReferenceSystems();
-        if (!componentsSrc.empty() &&
-            componentsSrc.size() == componentsDst.size()) {
-            if (componentsSrc[0]->extractGeographicCRS() &&
-                componentsDst[0]->extractGeographicCRS()) {
-
-                std::vector<CoordinateOperationNNPtr> verticalTransforms;
-                if (componentsSrc.size() >= 2 &&
-                    componentsSrc[1]->extractVerticalCRS() &&
-                    componentsDst[1]->extractVerticalCRS()) {
-                    verticalTransforms = createOperations(
-                        componentsSrc[1], componentsDst[1], context);
-                }
+// ---------------------------------------------------------------------------
 
-                for (const auto &verticalTransform : verticalTransforms) {
-                    auto interpolationGeogCRS =
-                        NN_NO_CHECK(componentsSrc[0]->extractGeographicCRS());
-                    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 =
-                                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<
-                                                crs::GeographicCRS>(
-                                        nn_interpTransformCRS));
-                            }
-                        }
-                    } else {
-                        auto compSrc0BoundCrs = dynamic_cast<crs::BoundCRS *>(
-                            componentsSrc[0].get());
-                        auto compDst0BoundCrs = dynamic_cast<crs::BoundCRS *>(
-                            componentsDst[0].get());
-                        if (compSrc0BoundCrs && compDst0BoundCrs &&
-                            dynamic_cast<crs::GeographicCRS *>(
-                                compSrc0BoundCrs->hubCRS().get()) &&
-                            compSrc0BoundCrs->hubCRS()->_isEquivalentTo(
-                                compDst0BoundCrs->hubCRS().get())) {
-                            interpolationGeogCRS =
-                                NN_NO_CHECK(util::nn_dynamic_pointer_cast<
-                                            crs::GeographicCRS>(
-                                    compSrc0BoundCrs->hubCRS()));
+void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::CompoundCRS *compoundSrc,
+    const crs::CompoundCRS *compoundDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    const auto &componentsSrc = compoundSrc->componentReferenceSystems();
+    const auto &componentsDst = compoundDst->componentReferenceSystems();
+    if (!componentsSrc.empty() &&
+        componentsSrc.size() == componentsDst.size()) {
+        if (componentsSrc[0]->extractGeographicCRS() &&
+            componentsDst[0]->extractGeographicCRS()) {
+
+            std::vector<CoordinateOperationNNPtr> verticalTransforms;
+            if (componentsSrc.size() >= 2 &&
+                componentsSrc[1]->extractVerticalCRS() &&
+                componentsDst[1]->extractVerticalCRS()) {
+                verticalTransforms = createOperations(
+                    componentsSrc[1], componentsDst[1], context);
+            }
+
+            for (const auto &verticalTransform : verticalTransforms) {
+                auto interpolationGeogCRS =
+                    NN_NO_CHECK(componentsSrc[0]->extractGeographicCRS());
+                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 = NN_NO_CHECK(
+                                util::nn_dynamic_pointer_cast<
+                                    crs::GeographicCRS>(nn_interpTransformCRS));
                         }
                     }
-                    auto opSrcCRSToGeogCRS = createOperations(
-                        componentsSrc[0], interpolationGeogCRS, context);
-                    auto opGeogCRStoDstCRS = createOperations(
-                        interpolationGeogCRS, componentsDst[0], context);
-                    for (const auto &opSrc : opSrcCRSToGeogCRS) {
-                        for (const auto &opDst : opGeogCRStoDstCRS) {
+                } else {
+                    auto compSrc0BoundCrs =
+                        dynamic_cast<crs::BoundCRS *>(componentsSrc[0].get());
+                    auto compDst0BoundCrs =
+                        dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
+                    if (compSrc0BoundCrs && compDst0BoundCrs &&
+                        dynamic_cast<crs::GeographicCRS *>(
+                            compSrc0BoundCrs->hubCRS().get()) &&
+                        compSrc0BoundCrs->hubCRS()->_isEquivalentTo(
+                            compDst0BoundCrs->hubCRS().get())) {
+                        interpolationGeogCRS = NN_NO_CHECK(
+                            util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
+                                compSrc0BoundCrs->hubCRS()));
+                    }
+                }
+                auto opSrcCRSToGeogCRS = createOperations(
+                    componentsSrc[0], interpolationGeogCRS, context);
+                auto opGeogCRStoDstCRS = createOperations(
+                    interpolationGeogCRS, componentsDst[0], context);
+                for (const auto &opSrc : opSrcCRSToGeogCRS) {
+                    for (const auto &opDst : opGeogCRStoDstCRS) {
 
-                            try {
-                                auto op = createHorizVerticalHorizPROJBased(
-                                    sourceCRS, targetCRS, opSrc,
-                                    verticalTransform, opDst,
-                                    interpolationGeogCRS, true);
-                                res.emplace_back(op);
-                            } catch (
-                                const InvalidOperationEmptyIntersection &) {
-                            }
+                        try {
+                            auto op = createHorizVerticalHorizPROJBased(
+                                sourceCRS, targetCRS, opSrc, verticalTransform,
+                                opDst, interpolationGeogCRS, true);
+                            res.emplace_back(op);
+                        } catch (const InvalidOperationEmptyIntersection &) {
                         }
                     }
                 }
+            }
 
-                if (verticalTransforms.empty()) {
-                    return createOperations(componentsSrc[0], componentsDst[0],
-                                            context);
-                }
+            if (verticalTransforms.empty()) {
+                res = createOperations(componentsSrc[0], componentsDst[0],
+                                       context);
             }
         }
     }
+}
 
-    // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to
-    // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx
-    // +type=crs'
-    if (boundSrc && compoundDst) {
-        const auto &componentsDst = compoundDst->componentReferenceSystems();
-        if (!componentsDst.empty()) {
-            auto compDst0BoundCrs =
-                dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
-            if (compDst0BoundCrs) {
-                auto boundSrcHubAsGeogCRS = dynamic_cast<crs::GeographicCRS *>(
-                    boundSrc->hubCRS().get());
-                auto compDst0BoundCrsHubAsGeogCRS =
-                    dynamic_cast<crs::GeographicCRS *>(
-                        compDst0BoundCrs->hubCRS().get());
-                if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) {
-                    const auto &boundSrcHubAsGeogCRSDatum =
-                        boundSrcHubAsGeogCRS->datum();
-                    const auto &compDst0BoundCrsHubAsGeogCRSDatum =
-                        compDst0BoundCrsHubAsGeogCRS->datum();
-                    if (boundSrcHubAsGeogCRSDatum &&
-                        compDst0BoundCrsHubAsGeogCRSDatum &&
-                        boundSrcHubAsGeogCRSDatum->_isEquivalentTo(
-                            compDst0BoundCrsHubAsGeogCRSDatum.get())) {
-                        auto cs = cs::EllipsoidalCS::
-                            createLatitudeLongitudeEllipsoidalHeight(
-                                common::UnitOfMeasure::DEGREE,
-                                common::UnitOfMeasure::METRE);
-                        auto intermGeog3DCRS = util::nn_static_pointer_cast<
-                            crs::CRS>(crs::GeographicCRS::create(
-                            util::PropertyMap()
-                                .set(common::IdentifiedObject::NAME_KEY,
-                                     boundSrcHubAsGeogCRS->nameStr())
-                                .set(
-                                    common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
-                                    metadata::Extent::WORLD),
-                            NN_NO_CHECK(boundSrcHubAsGeogCRSDatum), cs));
-                        auto sourceToGeog3DOps = createOperations(
-                            sourceCRS, intermGeog3DCRS, context);
-                        auto geog3DToTargetOps = createOperations(
-                            intermGeog3DCRS, targetCRS, context);
-                        for (const auto &opSrc : sourceToGeog3DOps) {
-                            for (const auto &opDst : geog3DToTargetOps) {
-                                if (opSrc->targetCRS() && opDst->sourceCRS() &&
-                                    !opSrc->targetCRS()->_isEquivalentTo(
-                                        opDst->sourceCRS().get())) {
-                                    // Shouldn't happen normally, but typically
-                                    // one of them can be 2D and the other 3D
-                                    // due to above createOperations() not
-                                    // exactly setting the expected source and
-                                    // target CRS.
-                                    // So create an adapter operation...
-                                    auto intermOps = createOperations(
-                                        NN_NO_CHECK(opSrc->targetCRS()),
-                                        NN_NO_CHECK(opDst->sourceCRS()),
-                                        context);
-                                    if (!intermOps.empty()) {
-                                        res.emplace_back(
-                                            ConcatenatedOperation::
-                                                createComputeMetadata(
-                                                    {opSrc, intermOps.front(),
-                                                     opDst},
-                                                    !allowEmptyIntersection));
-                                    }
-                                } else {
+// ---------------------------------------------------------------------------
+
+void CoordinateOperationFactory::Private::createOperationsBoundToCompound(
+    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    Private::Context &context, const crs::BoundCRS *boundSrc,
+    const crs::CompoundCRS *compoundDst,
+    std::vector<CoordinateOperationNNPtr> &res) {
+
+    const auto &componentsDst = compoundDst->componentReferenceSystems();
+    if (!componentsDst.empty()) {
+        auto compDst0BoundCrs =
+            dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
+        if (compDst0BoundCrs) {
+            auto boundSrcHubAsGeogCRS =
+                dynamic_cast<crs::GeographicCRS *>(boundSrc->hubCRS().get());
+            auto compDst0BoundCrsHubAsGeogCRS =
+                dynamic_cast<crs::GeographicCRS *>(
+                    compDst0BoundCrs->hubCRS().get());
+            if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) {
+                const auto &boundSrcHubAsGeogCRSDatum =
+                    boundSrcHubAsGeogCRS->datum();
+                const auto &compDst0BoundCrsHubAsGeogCRSDatum =
+                    compDst0BoundCrsHubAsGeogCRS->datum();
+                if (boundSrcHubAsGeogCRSDatum &&
+                    compDst0BoundCrsHubAsGeogCRSDatum &&
+                    boundSrcHubAsGeogCRSDatum->_isEquivalentTo(
+                        compDst0BoundCrsHubAsGeogCRSDatum.get())) {
+                    auto cs = cs::EllipsoidalCS::
+                        createLatitudeLongitudeEllipsoidalHeight(
+                            common::UnitOfMeasure::DEGREE,
+                            common::UnitOfMeasure::METRE);
+                    auto intermGeog3DCRS = util::nn_static_pointer_cast<
+                        crs::CRS>(crs::GeographicCRS::create(
+                        util::PropertyMap()
+                            .set(common::IdentifiedObject::NAME_KEY,
+                                 boundSrcHubAsGeogCRS->nameStr())
+                            .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                                 metadata::Extent::WORLD),
+                        NN_NO_CHECK(boundSrcHubAsGeogCRSDatum), cs));
+                    auto sourceToGeog3DOps =
+                        createOperations(sourceCRS, intermGeog3DCRS, context);
+                    auto geog3DToTargetOps =
+                        createOperations(intermGeog3DCRS, targetCRS, context);
+                    for (const auto &opSrc : sourceToGeog3DOps) {
+                        for (const auto &opDst : geog3DToTargetOps) {
+                            if (opSrc->targetCRS() && opDst->sourceCRS() &&
+                                !opSrc->targetCRS()->_isEquivalentTo(
+                                    opDst->sourceCRS().get())) {
+                                // Shouldn't happen normally, but typically
+                                // one of them can be 2D and the other 3D
+                                // due to above createOperations() not
+                                // exactly setting the expected source and
+                                // target CRS.
+                                // So create an adapter operation...
+                                auto intermOps = createOperations(
+                                    NN_NO_CHECK(opSrc->targetCRS()),
+                                    NN_NO_CHECK(opDst->sourceCRS()), context);
+                                if (!intermOps.empty()) {
                                     res.emplace_back(
                                         ConcatenatedOperation::
                                             createComputeMetadata(
-                                                {opSrc, opDst},
+                                                {opSrc, intermOps.front(),
+                                                 opDst},
                                                 !allowEmptyIntersection));
                                 }
+                            } else {
+                                res.emplace_back(
+                                    ConcatenatedOperation::
+                                        createComputeMetadata(
+                                            {opSrc, opDst},
+                                            !allowEmptyIntersection));
                             }
                         }
                     }
@@ -13617,13 +14044,6 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
     }
-
-    // reverse of previous case
-    if (boundDst && compoundSrc) {
-        return applyInverse(createOperations(targetCRS, sourceCRS, context));
-    }
-
-    return res;
 }
 //! @endcond
 
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 5d3af369..f16ed7cf 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -2170,7 +2170,7 @@ GeographicCRS::demoteTo2D(const std::string &newName,
     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
+        // First check if there is a Geographic 2D CRS in the database
         // of the same name.
         // This is the common practice in the EPSG dataset.
         if (dbContext && l_identifiers.size() == 1) {
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 993b4eeb..44e44dbe 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -2295,6 +2295,20 @@ AuthorityFactory::createConversion(const std::string &code) const {
 
     auto res = d->runWithCodeParam(sql, code);
     if (res.empty()) {
+        try {
+            // Conversions using methods Change of Vertical Unit or
+            // Height Depth Reveral are stored in other_transformation
+            auto op = createCoordinateOperation(
+                code, false /* allowConcatenated */,
+                false /* usePROJAlternativeGridNames */,
+                "other_transformation");
+            auto conv =
+                util::nn_dynamic_pointer_cast<operation::Conversion>(op);
+            if (conv) {
+                return NN_NO_CHECK(conv);
+            }
+        } catch (const std::exception &) {
+        }
         throw NoSuchAuthorityCodeException("conversion not found",
                                            d->authority(), code);
     }
@@ -3118,13 +3132,16 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
                     .set(metadata::Identifier::CODE_KEY, method_code)
                     .set(common::IdentifiedObject::NAME_KEY, method_name);
 
-            if (method_auth_name == metadata::Identifier::EPSG &&
-                operation::isAxisOrderReversal(
-                    std::atoi(method_code.c_str()))) {
-                auto op = operation::Conversion::create(props, propsMethod,
-                                                        parameters, values);
-                op->setCRSs(sourceCRS, targetCRS, nullptr);
-                return op;
+            if (method_auth_name == metadata::Identifier::EPSG) {
+                int method_code_int = std::atoi(method_code.c_str());
+                if (operation::isAxisOrderReversal(method_code_int) ||
+                    method_code_int == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT ||
+                    method_code_int == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) {
+                    auto op = operation::Conversion::create(props, propsMethod,
+                                                            parameters, values);
+                    op->setCRSs(sourceCRS, targetCRS, nullptr);
+                    return op;
+                }
             }
             return operation::Transformation::create(
                 props, sourceCRS, targetCRS, nullptr, propsMethod, parameters,
@@ -4603,6 +4620,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/src/iso19111/io.cpp b/src/iso19111/io.cpp
index fc66b6c9..50ad5a87 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3476,6 +3476,14 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
     }
     propertiesMethod.set(IdentifiedObject::NAME_KEY, projectionName);
 
+    std::vector<bool> foundParameters;
+    if (mapping) {
+        size_t countParams = 0;
+        while (mapping->params[countParams] != nullptr) {
+            ++countParams;
+        }
+        foundParameters.resize(countParams);
+    }
     for (const auto &childNode : projCRSNode->GP()->children()) {
         if (ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER)) {
             const auto &childNodeChildren = childNode->GP()->children();
@@ -3496,11 +3504,18 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
                     continue;
                 }
             }
-            const auto *paramMapping =
+            auto *paramMapping =
                 mapping ? getMappingFromWKT1(mapping, parameterName) : nullptr;
             if (mapping &&
                 mapping->epsg_code == EPSG_CODE_METHOD_MERCATOR_VARIANT_B &&
                 ci_equal(parameterName, "latitude_of_origin")) {
+                for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) {
+                    if (mapping->params[idx]->epsg_code ==
+                        EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN) {
+                        foundParameters[idx] = true;
+                        break;
+                    }
+                }
                 parameterName = EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN;
                 propertiesParameter.set(
                     Identifier::CODE_KEY,
@@ -3508,6 +3523,12 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
                 propertiesParameter.set(Identifier::CODESPACE_KEY,
                                         Identifier::EPSG);
             } else if (paramMapping) {
+                for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) {
+                    if (mapping->params[idx] == paramMapping) {
+                        foundParameters[idx] = true;
+                        break;
+                    }
+                }
                 parameterName = paramMapping->wkt2_name;
                 if (paramMapping->epsg_code != 0) {
                     propertiesParameter.set(Identifier::CODE_KEY,
@@ -3531,6 +3552,38 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
         }
     }
 
+    // Add back important parameters that should normally be present, but
+    // are sometimes missing. Currently we only deal with Scale factor at
+    // natural origin. This is to avoid a default value of 0 to slip in later.
+    // But such WKT should be considered invalid.
+    if (mapping) {
+        for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) {
+            if (!foundParameters[idx] &&
+                mapping->params[idx]->epsg_code ==
+                    EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) {
+
+                emitRecoverableWarning(
+                    "The WKT string lacks a value "
+                    "for " EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN
+                    ". Default it to 1.");
+
+                PropertyMap propertiesParameter;
+                propertiesParameter.set(
+                    Identifier::CODE_KEY,
+                    EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN);
+                propertiesParameter.set(Identifier::CODESPACE_KEY,
+                                        Identifier::EPSG);
+                propertiesParameter.set(
+                    IdentifiedObject::NAME_KEY,
+                    EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN);
+                parameters.push_back(
+                    OperationParameter::create(propertiesParameter));
+                values.push_back(ParameterValue::create(
+                    Measure(1.0, UnitOfMeasure::SCALE_UNITY)));
+            }
+        }
+    }
+
     return Conversion::create(
                PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"),
                propertiesMethod, parameters, values)
diff --git a/src/proj_constants.h b/src/proj_constants.h
index 619d9d91..62cf94be 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -592,4 +592,9 @@
 #define EPSG_NAME_METHOD_AXIS_ORDER_REVERSAL_3D \
     "Axis Order Reversal (Geographic3D horizontal)"
 
+/* ------------------------------------------------------------------------ */
+
+#define EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL 1068
+#define EPSG_NAME_METHOD_HEIGHT_DEPTH_REVERSAL "Height Depth Reversal"
+
 #endif /* PROJ_CONSTANTS_INCLUDED */