proj_crs_create_bound_crs_to_WGS84(): make it work on verticalCRS/compoundCRS such as EPSG:4326+5773 and EPSG:4326+3855

2 files changed, 110 insertions, 24 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 6593d3e9..6fcf4d30 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8702,15 +8702,6 @@ _getHeightToGeographic3DFilename(const Transformation *op, bool allowInverse) {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
-const std::string &Transformation::getHeightToGeographic3DFilename() const {
-
-    return _getHeightToGeographic3DFilename(this, false);
-}
-//! @endcond
-
-// ---------------------------------------------------------------------------
-
-//! @cond Doxygen_Suppress
 static bool
 isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
                                      bool allowInverse) {
@@ -8765,6 +8756,27 @@ isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+const std::string &Transformation::getHeightToGeographic3DFilename() const {
+
+    const std::string &ret = _getHeightToGeographic3DFilename(this, false);
+    if (!ret.empty())
+        return ret;
+    if (isGeographic3DToGravityRelatedHeight(method(), false)) {
+        const auto &fileParameter =
+            parameterValue(EPSG_NAME_PARAMETER_GEOID_CORRECTION_FILENAME,
+                           EPSG_CODE_PARAMETER_GEOID_CORRECTION_FILENAME);
+        if (fileParameter &&
+            fileParameter->type() == ParameterValue::Type::FILENAME) {
+            return fileParameter->valueFile();
+        }
+    }
+    return nullString;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 static util::PropertyMap
 createSimilarPropertiesMethod(common::IdentifiedObjectNNPtr obj) {
     util::PropertyMap map;
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index e96b3cc9..ecbd39e1 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -401,23 +401,20 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
         }
     }
 
-    auto geodCRS = util::nn_dynamic_pointer_cast<GeodeticCRS>(thisAsCRS);
-    auto geogCRS = extractGeographicCRS();
-    auto hubCRS = util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326);
-    if (geodCRS && !geogCRS) {
-        if (geodCRS->_isEquivalentTo(GeographicCRS::EPSG_4978.get(),
-                                     util::IComparable::Criterion::EQUIVALENT,
-                                     dbContext)) {
-            return thisAsCRS;
+    auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
+    if (compoundCRS) {
+        const auto &comps = compoundCRS->componentReferenceSystems();
+        if (comps.size() == 2) {
+            auto horiz = comps[0]->createBoundCRSToWGS84IfPossible(
+                dbContext, allowIntermediateCRSUse);
+            auto vert = comps[1]->createBoundCRSToWGS84IfPossible(
+                dbContext, allowIntermediateCRSUse);
+            if (horiz.get() != comps[0].get() || vert.get() != comps[1].get()) {
+                return CompoundCRS::create(createPropertyMap(this),
+                                           {horiz, vert});
+            }
         }
-        hubCRS = util::nn_static_pointer_cast<CRS>(GeodeticCRS::EPSG_4978);
-    } else if (!geogCRS ||
-               geogCRS->_isEquivalentTo(
-                   GeographicCRS::EPSG_4326.get(),
-                   util::IComparable::Criterion::EQUIVALENT, dbContext)) {
         return thisAsCRS;
-    } else {
-        geodCRS = geogCRS;
     }
 
     if (!dbContext) {
@@ -443,6 +440,83 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
     if (authorities.empty()) {
         authorities.emplace_back();
     }
+
+    // Vertical CRS ?
+    auto vertCRS = dynamic_cast<const VerticalCRS *>(this);
+    if (vertCRS) {
+        auto hubCRS =
+            util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4979);
+        for (const auto &authority : authorities) {
+            try {
+
+                auto authFactory = io::AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext),
+                    authority == "any" ? std::string() : authority);
+                auto ctxt = operation::CoordinateOperationContext::create(
+                    authFactory, extent, 0.0);
+                ctxt->setAllowUseIntermediateCRS(allowIntermediateCRSUse);
+                // ctxt->setSpatialCriterion(
+                //    operation::CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+                auto list = operation::CoordinateOperationFactory::create()
+                                ->createOperations(hubCRS, thisAsCRS, ctxt);
+                CRSPtr candidateBoundCRS;
+                for (const auto &op : list) {
+                    auto transf = util::nn_dynamic_pointer_cast<
+                        operation::Transformation>(op);
+                    // Only keep transformations that use a known grid
+                    if (transf && !transf->hasBallparkTransformation()) {
+                        auto gridsNeeded = transf->gridsNeeded(dbContext, true);
+                        bool gridsKnown = !gridsNeeded.empty();
+                        for (const auto &gridDesc : gridsNeeded) {
+                            if (gridDesc.packageName.empty() &&
+                                !(!gridDesc.url.empty() &&
+                                  gridDesc.openLicense) &&
+                                !gridDesc.available) {
+                                gridsKnown = false;
+                                break;
+                            }
+                        }
+                        if (gridsKnown) {
+                            if (candidateBoundCRS) {
+                                candidateBoundCRS = nullptr;
+                                break;
+                            }
+                            candidateBoundCRS =
+                                BoundCRS::create(thisAsCRS, hubCRS,
+                                                 NN_NO_CHECK(transf))
+                                    .as_nullable();
+                        }
+                    }
+                }
+                if (candidateBoundCRS) {
+                    return NN_NO_CHECK(candidateBoundCRS);
+                }
+            } catch (const std::exception &) {
+            }
+        }
+        return thisAsCRS;
+    }
+
+    // Geodetic/geographic CRS ?
+    auto geodCRS = util::nn_dynamic_pointer_cast<GeodeticCRS>(thisAsCRS);
+    auto geogCRS = extractGeographicCRS();
+    auto hubCRS = util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326);
+    if (geodCRS && !geogCRS) {
+        if (geodCRS->_isEquivalentTo(GeographicCRS::EPSG_4978.get(),
+                                     util::IComparable::Criterion::EQUIVALENT,
+                                     dbContext)) {
+            return thisAsCRS;
+        }
+        hubCRS = util::nn_static_pointer_cast<CRS>(GeodeticCRS::EPSG_4978);
+    } else if (!geogCRS ||
+               geogCRS->_isEquivalentTo(
+                   GeographicCRS::EPSG_4326.get(),
+                   util::IComparable::Criterion::EQUIVALENT, dbContext)) {
+        return thisAsCRS;
+    } else {
+        geodCRS = geogCRS;
+    }
+
     for (const auto &authority : authorities) {
         try {