[Backport 6.2] Improve vertical transforms (#1613)

[Backport 6.2] Improve vertical transforms

2 files changed, 361 insertions, 133 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 9d20b396..f0ebaebd 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,8 @@ struct CoordinateOperationFactory::Private {
         const crs::CRSNNPtr &targetCRS;
         const CoordinateOperationContextNNPtr &context;
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
+        bool inCreateOperationsGeogToVertWithIntermediate = false;
+        bool skipHorizontalTransformation = false;
 
         Context(const crs::CRSNNPtr &sourceCRSIn,
                 const crs::CRSNNPtr &targetCRSIn,
@@ -10193,6 +10197,11 @@ 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);
@@ -10309,16 +10318,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 +10685,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 +10950,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 +10990,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 +11032,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 +11045,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 +11067,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 +11156,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 +11169,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 +12077,51 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
 // ---------------------------------------------------------------------------
 
+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) {
@@ -12106,11 +12238,16 @@ CoordinateOperationFactory::Private::createOperations(
 
     auto geodSrc = dynamic_cast<const crs::GeodeticCRS *>(sourceCRS.get());
     auto geodDst = dynamic_cast<const crs::GeodeticCRS *>(targetCRS.get());
+    auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+    auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
+    auto vertSrc = dynamic_cast<const crs::VerticalCRS *>(sourceCRS.get());
+    auto vertDst = dynamic_cast<const crs::VerticalCRS *>(targetCRS.get());
 
     // First look-up if the registry provide us with operations.
     auto derivedSrc = dynamic_cast<const crs::DerivedCRS *>(sourceCRS.get());
     auto derivedDst = dynamic_cast<const crs::DerivedCRS *>(targetCRS.get());
-    if (context.context->getAuthorityFactory() &&
+    auto authFactory = context.context->getAuthorityFactory();
+    if (authFactory &&
         (derivedSrc == nullptr ||
          !derivedSrc->baseCRS()->_isEquivalentTo(
              targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) &&
@@ -12144,6 +12281,21 @@ CoordinateOperationFactory::Private::createOperations(
                         util::IComparable::Criterion::EQUIVALENT);
             }
 
+            // 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) {
@@ -12202,10 +12354,6 @@ CoordinateOperationFactory::Private::createOperations(
                 "celestial body");
         }
 
-        auto geogSrc =
-            dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
-        auto geogDst =
-            dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
         if (geogSrc && geogDst) {
             return createOperationsGeogToGeog(res, sourceCRS, targetCRS,
                                               geogSrc, geogDst);
@@ -12293,7 +12441,6 @@ CoordinateOperationFactory::Private::createOperations(
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
-    auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
     if (boundSrc && geogDst) {
         const auto &hubSrc = boundSrc->hubCRS();
         auto hubSrcGeog =
@@ -12453,27 +12600,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;
+                    }
                 }
             }
         }
@@ -12482,14 +12634,12 @@ CoordinateOperationFactory::Private::createOperations(
     }
 
     // reverse of previous case
-    auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
     if (geogSrc && boundDst) {
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
     // vertCRS (as boundCRS with transformation to target vertCRS) to
     // vertCRS
-    auto vertDst = dynamic_cast<const crs::VerticalCRS *>(targetCRS.get());
     if (boundSrc && vertDst) {
         auto baseSrcVert =
             dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
@@ -12506,7 +12656,6 @@ CoordinateOperationFactory::Private::createOperations(
     }
 
     // reverse of previous case
-    auto vertSrc = dynamic_cast<const crs::VerticalCRS *>(sourceCRS.get());
     if (boundDst && vertSrc) {
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
@@ -12550,7 +12699,6 @@ CoordinateOperationFactory::Private::createOperations(
     if (vertSrc && geogDst) {
 
         if (vertSrc->identifiers().empty()) {
-            const auto &authFactory = context.context->getAuthorityFactory();
             const auto &vertSrcName = vertSrc->nameStr();
             if (authFactory != nullptr && vertSrcName != "unnamed" &&
                 vertSrcName != "unknown") {
@@ -12699,8 +12847,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) {
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 1f40f1f0..0f6790c0 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -3312,9 +3312,9 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes(
     bool discardSuperseded) const {
 
     auto cacheKey(d->authority());
-    cacheKey += sourceCRSAuthName;
+    cacheKey += sourceCRSAuthName.empty() ? "{empty}" : sourceCRSAuthName;
     cacheKey += sourceCRSCode;
-    cacheKey += targetCRSAuthName;
+    cacheKey += targetCRSAuthName.empty() ? "{empty}" : targetCRSAuthName;
     cacheKey += targetCRSCode;
     cacheKey += (usePROJAlternativeGridNames ? '1' : '0');
     cacheKey += (discardIfMissingGrid ? '1' : '0');
@@ -3326,27 +3326,30 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes(
         return list;
     }
 
-    // Look-up first for conversion which is the most precise.
-    std::string sql("SELECT conversion_auth_name, "
-                    "geodetic_crs_auth_name, geodetic_crs_code FROM "
-                    "projected_crs WHERE auth_name = ? AND code = ?");
-    auto params = ListOfParams{targetCRSAuthName, targetCRSCode};
-    auto res = d->run(sql, params);
-    if (!res.empty()) {
-        const auto &row = res.front();
-        bool ok = row[1] == sourceCRSAuthName && row[2] == sourceCRSCode;
-        if (ok && d->hasAuthorityRestriction()) {
-            ok = row[0] == d->authority();
-        }
-        if (ok) {
-            auto targetCRS = d->createFactory(targetCRSAuthName)
-                                 ->createProjectedCRS(targetCRSCode);
-            auto conv = targetCRS->derivingConversion();
-            list.emplace_back(conv);
-            d->context()->d->cache(cacheKey, list);
-            return list;
+    if (!targetCRSAuthName.empty()) {
+        // Look-up first for conversion which is the most precise.
+        std::string sql("SELECT conversion_auth_name, "
+                        "geodetic_crs_auth_name, geodetic_crs_code FROM "
+                        "projected_crs WHERE auth_name = ? AND code = ?");
+        auto params = ListOfParams{targetCRSAuthName, targetCRSCode};
+        auto res = d->run(sql, params);
+        if (!res.empty()) {
+            const auto &row = res.front();
+            bool ok = row[1] == sourceCRSAuthName && row[2] == sourceCRSCode;
+            if (ok && d->hasAuthorityRestriction()) {
+                ok = row[0] == d->authority();
+            }
+            if (ok) {
+                auto targetCRS = d->createFactory(targetCRSAuthName)
+                                     ->createProjectedCRS(targetCRSCode);
+                auto conv = targetCRS->derivingConversion();
+                list.emplace_back(conv);
+                d->context()->d->cache(cacheKey, list);
+                return list;
+            }
         }
     }
+    std::string sql;
     if (discardSuperseded) {
         sql = "SELECT cov.auth_name, cov.code, cov.table_name, "
               "ss.replacement_auth_name, ss.replacement_code FROM "
@@ -3358,20 +3361,26 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes(
               "ss.superseded_auth_name = cov.auth_name AND "
               "ss.superseded_code = cov.code AND "
               "ss.superseded_table_name = ss.replacement_table_name "
-              "WHERE source_crs_auth_name = ? AND source_crs_code = ? AND "
-              "target_crs_auth_name = ? AND target_crs_code = ? AND "
-              "cov.deprecated = 0";
+              "WHERE ";
     } else {
         sql = "SELECT cov.auth_name, cov.code, cov.table_name FROM "
               "coordinate_operation_view cov JOIN area "
               "ON cov.area_of_use_auth_name = area.auth_name AND "
               "cov.area_of_use_code = area.code "
-              "WHERE source_crs_auth_name = ? AND source_crs_code = ? AND "
-              "target_crs_auth_name = ? AND target_crs_code = ? AND "
-              "cov.deprecated = 0";
+              "WHERE ";
+    }
+    ListOfParams params;
+    if (!sourceCRSAuthName.empty()) {
+        sql += "source_crs_auth_name = ? AND source_crs_code = ? AND ";
+        params.emplace_back(sourceCRSAuthName);
+        params.emplace_back(sourceCRSCode);
+    }
+    if (!targetCRSAuthName.empty()) {
+        sql += "target_crs_auth_name = ? AND target_crs_code = ? AND ";
+        params.emplace_back(targetCRSAuthName);
+        params.emplace_back(targetCRSCode);
     }
-    params = {sourceCRSAuthName, sourceCRSCode, targetCRSAuthName,
-              targetCRSCode};
+    sql += "cov.deprecated = 0";
     if (d->hasAuthorityRestriction()) {
         sql += " AND cov.auth_name = ?";
         params.emplace_back(d->authority());
@@ -3379,7 +3388,7 @@ AuthorityFactory::createFromCoordinateReferenceSystemCodes(
     sql += " ORDER BY pseudo_area_from_swne(south_lat, west_lon, north_lat, "
            "east_lon) DESC, "
            "(CASE WHEN accuracy is NULL THEN 1 ELSE 0 END), accuracy";
-    res = d->run(sql, params);
+    auto res = d->run(sql, params);
     std::set<std::pair<std::string, std::string>> setTransf;
     if (discardSuperseded) {
         for (const auto &row : res) {