Merge pull request #1608 from rouault/improve_transforms_fromto_wgs84_gXXXX_realizations

Improvements in transformations from/to WGS 84 (Gxxxx) realizations and vertical <--> geog transormations

1 files changed, 450 insertions, 94 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 768ef76f..aea8400c 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -56,7 +56,9 @@
 #include <string>
 #include <vector>
 
-#ifdef DEBUG
+// #define DEBUG
+// #define DEBUG_SORT
+#if defined(DEBUG) || defined(DEBUG_SORT)
 #include <iostream>
 #endif
 
@@ -10169,6 +10171,9 @@ struct CoordinateOperationFactory::Private {
         const crs::CRSNNPtr &targetCRS;
         const CoordinateOperationContextNNPtr &context;
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
+        bool inCreateOperationsThroughPreferredHub = false;
+        bool inCreateOperationsGeogToVertWithIntermediate = false;
+        bool skipHorizontalTransformation = false;
 
         Context(const crs::CRSNNPtr &sourceCRSIn,
                 const crs::CRSNNPtr &targetCRSIn,
@@ -10193,6 +10198,17 @@ struct CoordinateOperationFactory::Private {
         const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst,
         Context &context);
 
+    static void createOperationsThroughPreferredHub(
+        std::vector<CoordinateOperationNNPtr> &res,
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        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);
@@ -10309,16 +10325,12 @@ struct SortFunction {
             return false;
         }
 
-        if (iterA->second.hasGrids_ && iterB->second.hasGrids_) {
-            // Operations where grids are all available go before other
-            if (iterA->second.gridsAvailable_ &&
-                !iterB->second.gridsAvailable_) {
-                return true;
-            }
-            if (iterB->second.gridsAvailable_ &&
-                !iterA->second.gridsAvailable_) {
-                return false;
-            }
+        // Operations where grids are all available go before other
+        if (iterA->second.gridsAvailable_ && !iterB->second.gridsAvailable_) {
+            return true;
+        }
+        if (iterB->second.gridsAvailable_ && !iterA->second.gridsAvailable_) {
+            return false;
         }
 
         // Operations where grids are all known in our DB go before other
@@ -10680,7 +10692,35 @@ struct FilterResults {
         }
 
         // Sort !
-        std::sort(res.begin(), res.end(), SortFunction(map));
+        SortFunction sortFunc(map);
+        std::sort(res.begin(), res.end(), sortFunc);
+
+// Debug code to check consistency of the sort function
+#ifdef DEBUG_SORT
+        constexpr bool debugSort = true;
+#elif !defined(NDEBUG)
+        const bool debugSort = getenv("PROJ_DEBUG_SORT_FUNCT") != nullptr;
+#endif
+#if defined(DEBUG_SORT) || !defined(NDEBUG)
+        if (debugSort) {
+            const bool assertIfIssue =
+                !(getenv("PROJ_DEBUG_SORT_FUNCT_ASSERT") != nullptr);
+            for (size_t i = 0; i < res.size(); ++i) {
+                for (size_t j = i + 1; j < res.size(); ++j) {
+                    if (sortFunc(res[j], res[i])) {
+#ifdef DEBUG_SORT
+                        std::cerr << "Sorting issue with entry " << i << "("
+                                  << res[i]->nameStr() << ") and " << j << "("
+                                  << res[j]->nameStr() << ")" << std::endl;
+#endif
+                        if (assertIfIssue) {
+                            assert(false);
+                        }
+                    }
+                }
+            }
+        }
+#endif
     }
 
     // ----------------------------------------------------------------------
@@ -10917,18 +10957,21 @@ applyInverse(const std::vector<CoordinateOperationNNPtr> &list) {
 
 //! @cond Doxygen_Suppress
 
-static void buildSourceAndTargetCRSIds(
-    const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-    const CoordinateOperationContextNNPtr &context,
-    std::list<std::pair<std::string, std::string>> &sourceIds,
-    std::list<std::pair<std::string, std::string>> &targetIds) {
+static void buildCRSIds(const crs::CRSNNPtr &crs,
+                        const CoordinateOperationContextNNPtr &context,
+                        std::list<std::pair<std::string, std::string>> &ids) {
     const auto &authFactory = context->getAuthorityFactory();
     assert(authFactory);
     const auto &authFactoryName = authFactory->getAuthority();
 
-    const auto findObjectInDB = [&authFactory, &authFactoryName](
-        const crs::CRSNNPtr &crs,
-        std::list<std::pair<std::string, std::string>> &idList) {
+    for (const auto &id : crs->identifiers()) {
+        const auto &authName = *(id->codeSpace());
+        const auto &code = id->code();
+        if (!authName.empty()) {
+            ids.emplace_back(authName, code);
+        }
+    }
+    if (ids.empty()) {
         try {
             const auto tmpAuthFactory = io::AuthorityFactory::create(
                 authFactory->databaseContext(),
@@ -10954,35 +10997,37 @@ static void buildSourceAndTargetCRSIds(
                         matches.front().get(),
                         util::IComparable::Criterion::EQUIVALENT) &&
                     !matches.front()->identifiers().empty()) {
-                    const auto &ids = matches.front()->identifiers();
-                    idList.emplace_back(*(ids[0]->codeSpace()), ids[0]->code());
+                    const auto &tmpIds = matches.front()->identifiers();
+                    ids.emplace_back(*(tmpIds[0]->codeSpace()),
+                                     tmpIds[0]->code());
                 }
             }
         } catch (const std::exception &) {
         }
-    };
-
-    for (const auto &id : sourceCRS->identifiers()) {
-        const auto &authName = *(id->codeSpace());
-        const auto &code = id->code();
-        if (!authName.empty()) {
-            sourceIds.emplace_back(authName, code);
-        }
-    }
-    if (sourceIds.empty()) {
-        findObjectInDB(sourceCRS, sourceIds);
     }
+}
 
-    for (const auto &id : targetCRS->identifiers()) {
-        const auto &authName = *(id->codeSpace());
-        const auto &code = id->code();
-        if (!authName.empty()) {
-            targetIds.emplace_back(authName, code);
-        }
+// ---------------------------------------------------------------------------
+
+static std::vector<std::string>
+getCandidateAuthorities(const io::AuthorityFactoryPtr &authFactory,
+                        const std::string &srcAuthName,
+                        const std::string &targetAuthName) {
+    const auto &authFactoryName = authFactory->getAuthority();
+    std::vector<std::string> authorities;
+    if (authFactoryName == "any") {
+        authorities.emplace_back();
     }
-    if (targetIds.empty()) {
-        findObjectInDB(targetCRS, targetIds);
+    if (authFactoryName.empty()) {
+        authorities = authFactory->databaseContext()->getAllowedAuthorities(
+            srcAuthName, targetAuthName);
+        if (authorities.empty()) {
+            authorities.emplace_back();
+        }
+    } else {
+        authorities.emplace_back(authFactoryName);
     }
+    return authorities;
 }
 
 // ---------------------------------------------------------------------------
@@ -10994,12 +11039,11 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS,
                         const CoordinateOperationContextNNPtr &context) {
     const auto &authFactory = context->getAuthorityFactory();
     assert(authFactory);
-    const auto &authFactoryName = authFactory->getAuthority();
 
     std::list<std::pair<std::string, std::string>> sourceIds;
     std::list<std::pair<std::string, std::string>> targetIds;
-    buildSourceAndTargetCRSIds(sourceCRS, targetCRS, context, sourceIds,
-                               targetIds);
+    buildCRSIds(sourceCRS, context, sourceIds);
+    buildCRSIds(targetCRS, context, targetIds);
 
     for (const auto &idSrc : sourceIds) {
         const auto &srcAuthName = idSrc.first;
@@ -11008,20 +11052,8 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS,
             const auto &targetAuthName = idTarget.first;
             const auto &targetCode = idTarget.second;
 
-            std::vector<std::string> authorities;
-            if (authFactoryName == "any") {
-                authorities.emplace_back();
-            }
-            if (authFactoryName.empty()) {
-                authorities =
-                    authFactory->databaseContext()->getAllowedAuthorities(
-                        srcAuthName, targetAuthName);
-                if (authorities.empty()) {
-                    authorities.emplace_back();
-                }
-            } else {
-                authorities.emplace_back(authFactoryName);
-            }
+            const auto authorities(getCandidateAuthorities(
+                authFactory, srcAuthName, targetAuthName));
             for (const auto &authority : authorities) {
                 const auto tmpAuthFactory = io::AuthorityFactory::create(
                     authFactory->databaseContext(),
@@ -11042,6 +11074,81 @@ findOpsInRegistryDirect(const crs::CRSNNPtr &sourceCRS,
     }
     return std::vector<CoordinateOperationNNPtr>();
 }
+
+// ---------------------------------------------------------------------------
+
+// Look in the authority registry for operations from sourceCRS
+static std::vector<CoordinateOperationNNPtr>
+findOpsInRegistryDirectFrom(const crs::CRSNNPtr &sourceCRS,
+                            const CoordinateOperationContextNNPtr &context) {
+    const auto &authFactory = context->getAuthorityFactory();
+    assert(authFactory);
+
+    std::list<std::pair<std::string, std::string>> ids;
+    buildCRSIds(sourceCRS, context, ids);
+
+    for (const auto &id : ids) {
+        const auto &srcAuthName = id.first;
+        const auto &srcCode = id.second;
+
+        const auto authorities(
+            getCandidateAuthorities(authFactory, srcAuthName, srcAuthName));
+        for (const auto &authority : authorities) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(),
+                authority == "any" ? std::string() : authority);
+            auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
+                srcAuthName, srcCode, std::string(), std::string(),
+                context->getUsePROJAlternativeGridNames(),
+                context->getGridAvailabilityUse() ==
+                    CoordinateOperationContext::GridAvailabilityUse::
+                        DISCARD_OPERATION_IF_MISSING_GRID,
+                context->getDiscardSuperseded());
+            if (!res.empty()) {
+                return res;
+            }
+        }
+    }
+    return std::vector<CoordinateOperationNNPtr>();
+}
+
+// ---------------------------------------------------------------------------
+
+// Look in the authority registry for operations to targetCRS
+static std::vector<CoordinateOperationNNPtr>
+findOpsInRegistryDirectTo(const crs::CRSNNPtr &targetCRS,
+                          const CoordinateOperationContextNNPtr &context) {
+    const auto &authFactory = context->getAuthorityFactory();
+    assert(authFactory);
+
+    std::list<std::pair<std::string, std::string>> ids;
+    buildCRSIds(targetCRS, context, ids);
+
+    for (const auto &id : ids) {
+        const auto &targetAuthName = id.first;
+        const auto &targetCode = id.second;
+
+        const auto authorities(getCandidateAuthorities(
+            authFactory, targetAuthName, targetAuthName));
+        for (const auto &authority : authorities) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(),
+                authority == "any" ? std::string() : authority);
+            auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
+                std::string(), std::string(), targetAuthName, targetCode,
+                context->getUsePROJAlternativeGridNames(),
+                context->getGridAvailabilityUse() ==
+                    CoordinateOperationContext::GridAvailabilityUse::
+                        DISCARD_OPERATION_IF_MISSING_GRID,
+                context->getDiscardSuperseded());
+            if (!res.empty()) {
+                return res;
+            }
+        }
+    }
+    return std::vector<CoordinateOperationNNPtr>();
+}
+
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -11056,12 +11163,11 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate(
 
     const auto &authFactory = context->getAuthorityFactory();
     assert(authFactory);
-    const auto &authFactoryName = authFactory->getAuthority();
 
     std::list<std::pair<std::string, std::string>> sourceIds;
     std::list<std::pair<std::string, std::string>> targetIds;
-    buildSourceAndTargetCRSIds(sourceCRS, targetCRS, context, sourceIds,
-                               targetIds);
+    buildCRSIds(sourceCRS, context, sourceIds);
+    buildCRSIds(targetCRS, context, targetIds);
 
     for (const auto &idSrc : sourceIds) {
         const auto &srcAuthName = idSrc.first;
@@ -11070,20 +11176,8 @@ static std::vector<CoordinateOperationNNPtr> findsOpsInRegistryWithIntermediate(
             const auto &targetAuthName = idTarget.first;
             const auto &targetCode = idTarget.second;
 
-            std::vector<std::string> authorities;
-            if (authFactoryName == "any") {
-                authorities.emplace_back();
-            }
-            if (authFactoryName.empty()) {
-                authorities =
-                    authFactory->databaseContext()->getAllowedAuthorities(
-                        srcAuthName, targetAuthName);
-                if (authorities.empty()) {
-                    authorities.emplace_back();
-                }
-            } else {
-                authorities.emplace_back(authFactoryName);
-            }
+            const auto authorities(getCandidateAuthorities(
+                authFactory, srcAuthName, targetAuthName));
             for (const auto &authority : authorities) {
                 const auto tmpAuthFactory = io::AuthorityFactory::create(
                     authFactory->databaseContext(),
@@ -11990,6 +12084,167 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
 // ---------------------------------------------------------------------------
 
+void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub(
+    std::vector<CoordinateOperationNNPtr> &res, const crs::CRSNNPtr &sourceCRS,
+    const crs::CRSNNPtr &targetCRS, const crs::GeodeticCRS *geodSrc,
+    const crs::GeodeticCRS *geodDst, Private::Context &context) {
+
+    const auto &srcDatum = geodSrc->datum();
+    const auto &dstDatum = geodDst->datum();
+
+    if (!srcDatum || !dstDatum)
+        return;
+    const auto &srcDatumIds = srcDatum->identifiers();
+    const auto &dstDatumIds = dstDatum->identifiers();
+    if (srcDatumIds.empty() || dstDatumIds.empty())
+        return;
+
+    const auto &authFactory = context.context->getAuthorityFactory();
+
+    const auto srcAuthFactory = io::AuthorityFactory::create(
+        authFactory->databaseContext(), *(srcDatumIds.front()->codeSpace()));
+    const auto srcPreferredHubs =
+        srcAuthFactory->getPreferredHubGeodeticReferenceFrames(
+            srcDatumIds.front()->code());
+
+    const auto dstAuthFactory = io::AuthorityFactory::create(
+        authFactory->databaseContext(), *(dstDatumIds.front()->codeSpace()));
+    const auto dstPreferredHubs =
+        dstAuthFactory->getPreferredHubGeodeticReferenceFrames(
+            dstDatumIds.front()->code());
+    if (srcPreferredHubs.empty() && dstPreferredHubs.empty())
+        return;
+
+    // Currently if we have prefered hubs for both source and target, we
+    // will use only the one for target, arbitrarily... We could use boths
+    // but that would complicate things.
+    if (!srcPreferredHubs.empty() && dstPreferredHubs.empty()) {
+        std::vector<CoordinateOperationNNPtr> resTmp;
+        createOperationsThroughPreferredHub(resTmp, targetCRS, sourceCRS,
+                                            geodDst, geodSrc, context);
+        if (!resTmp.empty()) {
+            resTmp = applyInverse(resTmp);
+            res.insert(res.end(), resTmp.begin(), resTmp.end());
+        }
+        return;
+    }
+    assert(!dstPreferredHubs.empty());
+
+#ifdef DEBUG
+    EnterDebugLevel enterFunction;
+    debugTrace("createOperationsThroughPreferredHub(" +
+               objectAsStr(sourceCRS.get()) + "," +
+               objectAsStr(targetCRS.get()) + ")");
+#endif
+
+    struct AntiRecursionGuard {
+        Context &context;
+
+        explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
+            assert(!context.inCreateOperationsThroughPreferredHub);
+            context.inCreateOperationsThroughPreferredHub = true;
+        }
+
+        ~AntiRecursionGuard() {
+            context.inCreateOperationsThroughPreferredHub = false;
+        }
+    };
+    AntiRecursionGuard guard(context);
+
+    std::vector<crs::CRSNNPtr> candidatesIntermCRS;
+    for (const auto &datumHub : dstPreferredHubs) {
+        auto candidates =
+            findCandidateGeodCRSForDatum(authFactory, datumHub.get());
+        bool addedGeog2D = false;
+        for (const auto &intermCRS : candidates) {
+            auto geogCRS = dynamic_cast<crs::GeographicCRS *>(intermCRS.get());
+            if (geogCRS &&
+                geogCRS->coordinateSystem()->axisList().size() == 2) {
+                candidatesIntermCRS.emplace_back(intermCRS);
+                addedGeog2D = true;
+                break;
+            }
+        }
+        if (!addedGeog2D) {
+            candidatesIntermCRS.insert(candidatesIntermCRS.end(),
+                                       candidates.begin(), candidates.end());
+        }
+    }
+
+    const bool allowEmptyIntersection = true;
+    for (const auto &intermCRS : candidatesIntermCRS) {
+#ifdef DEBUG
+        EnterDebugLevel loopLevel;
+        debugTrace("try " + objectAsStr(sourceCRS.get()) + "->" +
+                   objectAsStr(intermCRS.get()) + "->" +
+                   objectAsStr(targetCRS.get()) + ")");
+        EnterDebugLevel loopLevel2;
+#endif
+        const auto opsFirst = createOperations(sourceCRS, intermCRS, context);
+        const auto opsLast = createOperations(intermCRS, targetCRS, context);
+        for (const auto &opFirst : opsFirst) {
+            for (const auto &opLast : opsLast) {
+                if (!opFirst->hasBallparkTransformation() ||
+                    !opLast->hasBallparkTransformation()) {
+                    try {
+                        res.emplace_back(
+                            ConcatenatedOperation::createComputeMetadata(
+                                {opFirst, opLast}, !allowEmptyIntersection));
+                    } catch (const InvalidOperationEmptyIntersection &) {
+                    }
+                }
+            }
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+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) {
@@ -12201,6 +12456,21 @@ CoordinateOperationFactory::Private::createOperations(
                 }
             }
 
+            // 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));
+            }
+
             if (res.empty() && !sameGeodeticDatum &&
                 !context.inCreateOperationsWithDatumPivotAntiRecursion &&
                 geodSrc && geodDst) {
@@ -12236,7 +12506,17 @@ CoordinateOperationFactory::Private::createOperations(
                     sourceCRS, targetCRS, context.context);
                 res.insert(res.end(), resWithIntermediate.begin(),
                            resWithIntermediate.end());
-                doFilterAndCheckPerfectOp = true;
+                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();
+                }
             }
         }
 
@@ -12505,27 +12785,32 @@ CoordinateOperationFactory::Private::createOperations(
             hubSrcGeog->coordinateSystem()->axisList().size() == 3 &&
             geogDst->coordinateSystem()->axisList().size() == 3) {
             auto opsFirst = createOperations(sourceCRS, hubSrc, context);
-            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 (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;
+                    if (!res.empty()) {
+                        return res;
+                    }
                 }
             }
         }
@@ -12747,8 +13032,79 @@ CoordinateOperationFactory::Private::createOperations(
             std::vector<CoordinateOperationNNPtr> verticalTransforms;
             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, context);
+                bool foundRegisteredTransformWithAllGridsAvailable = false;
+                for (const auto &op : verticalTransforms) {
+                    if (!op->identifiers().empty() && authFactory) {
+                        bool missingGrid = false;
+                        const auto gridsNeeded =
+                            op->gridsNeeded(authFactory->databaseContext());
+                        for (const auto &gridDesc : gridsNeeded) {
+                            if (!gridDesc.available) {
+                                missingGrid = 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());
+                            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());
+                    }
+                }
             }
             if (!horizTransforms.empty() && !verticalTransforms.empty()) {
                 for (const auto &horizTransform : horizTransforms) {