Merge pull request #2229 from rouault/fix_2228

Allow importing WKT1 COMPD_CS with a VERT_DATUM[Ellipsoid,2002]

3 files changed, 135 insertions, 20 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 88db81f9..bad7deea 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -94,7 +94,11 @@ struct CRS::Private {
     BoundCRSPtr canonicalBoundCRS_{};
     std::string extensionProj4_{};
     bool implicitCS_ = false;
+
     bool allowNonConformantWKT1Export_ = false;
+    // for what was initially a COMPD_CS with a VERT_CS with a datum type ==
+    // ellipsoidal height / 2002
+    VerticalCRSPtr originalVertCRS_{};
 
     void setImplicitCS(const util::PropertyMap &properties) {
         const auto pVal = properties.get("IMPLICIT_CS");
@@ -562,7 +566,13 @@ CRSNNPtr CRS::shallowClone() const { return _shallowClone(); }
 //! @cond Doxygen_Suppress
 
 CRSNNPtr CRS::allowNonConformantWKT1Export() const {
-    auto crs = shallowClone();
+    const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
+    if (boundCRS) {
+        return BoundCRS::create(
+            boundCRS->baseCRS()->allowNonConformantWKT1Export(),
+            boundCRS->hubCRS(), boundCRS->transformation());
+    }
+    auto crs(shallowClone());
     crs->d->allowNonConformantWKT1Export_ = true;
     return crs;
 }
@@ -573,6 +583,26 @@ CRSNNPtr CRS::allowNonConformantWKT1Export() const {
 
 //! @cond Doxygen_Suppress
 
+CRSNNPtr CRS::attachOriginalVertCRS(const VerticalCRSNNPtr &vertCRS) const {
+
+    const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
+    if (boundCRS) {
+        return BoundCRS::create(
+            boundCRS->baseCRS()->attachOriginalVertCRS(vertCRS),
+            boundCRS->hubCRS(), boundCRS->transformation());
+    }
+
+    auto crs(shallowClone());
+    crs->d->originalVertCRS_ = vertCRS.as_nullable();
+    return crs;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
+
 CRSNNPtr CRS::alterName(const std::string &newName) const {
     auto crs = shallowClone();
     auto newNameMod(newName);
@@ -1381,15 +1411,42 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     if (!isWKT2 && formatter->isStrict() && isGeographic &&
         axisList.size() != 2 &&
         oldAxisOutputRule != io::WKTFormatter::OutputAxisRule::NO) {
+
+        auto geogCRS2D = demoteTo2D(std::string(), dbContext);
+        if (dbContext) {
+            const auto res = geogCRS2D->identify(
+                io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "EPSG"));
+            if (res.size() == 1) {
+                const auto &front = res.front();
+                if (front.second == 100) {
+                    geogCRS2D = front.first;
+                }
+            }
+        }
+
         if (CRS::getPrivate()->allowNonConformantWKT1Export_) {
             formatter->startNode(io::WKTConstants::COMPD_CS, false);
             formatter->addQuotedString(l_name + " + " + l_name);
-            auto geogCRS = demoteTo2D(std::string(), dbContext);
-            geogCRS->_exportToWKT(formatter);
-            geogCRS->_exportToWKT(formatter);
+            geogCRS2D->_exportToWKT(formatter);
+            const auto oldTOWGSParameters = formatter->getTOWGS84Parameters();
+            formatter->setTOWGS84Parameters({});
+            geogCRS2D->_exportToWKT(formatter);
+            formatter->setTOWGS84Parameters(oldTOWGSParameters);
             formatter->endNode();
             return;
         }
+
+        auto &originalVertCRS = CRS::getPrivate()->originalVertCRS_;
+        if (originalVertCRS) {
+            formatter->startNode(io::WKTConstants::COMPD_CS, false);
+            formatter->addQuotedString(l_name + " + " +
+                                       originalVertCRS->nameStr());
+            geogCRS2D->_exportToWKT(formatter);
+            originalVertCRS->_exportToWKT(formatter);
+            formatter->endNode();
+            return;
+        }
+
         io::FormattingException::Throw(
             "WKT1 does not support Geographic 3D CRS.");
     }
@@ -2047,6 +2104,7 @@ GeodeticCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
 //! @cond Doxygen_Suppress
 struct GeographicCRS::Private {
     cs::EllipsoidalCSNNPtr coordinateSystem_;
+
     explicit Private(const cs::EllipsoidalCSNNPtr &csIn)
         : coordinateSystem_(csIn) {}
 };
@@ -3210,22 +3268,22 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     const auto &l_coordinateSystem = d->coordinateSystem();
     const auto &axisList = l_coordinateSystem->axisList();
     if (axisList.size() == 3 && !(isWKT2 && formatter->use2019Keywords())) {
+        auto projCRS2D = demoteTo2D(std::string(), dbContext);
+        if (dbContext) {
+            const auto res = projCRS2D->identify(
+                io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "EPSG"));
+            if (res.size() == 1) {
+                const auto &front = res.front();
+                if (front.second == 100) {
+                    projCRS2D = front.first;
+                }
+            }
+        }
+
         if (!formatter->useESRIDialect() &&
             CRS::getPrivate()->allowNonConformantWKT1Export_) {
             formatter->startNode(io::WKTConstants::COMPD_CS, false);
             formatter->addQuotedString(l_name + " + " + baseCRS()->nameStr());
-            auto projCRS2D = demoteTo2D(std::string(), dbContext);
-            if (dbContext) {
-                const auto res =
-                    projCRS2D->identify(io::AuthorityFactory::create(
-                        NN_NO_CHECK(dbContext), "EPSG"));
-                if (res.size() == 1) {
-                    const auto &front = res.front();
-                    if (front.second == 100) {
-                        projCRS2D = front.first;
-                    }
-                }
-            }
             projCRS2D->_exportToWKT(formatter);
             baseCRS()
                 ->demoteTo2D(std::string(), dbContext)
@@ -3233,6 +3291,18 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             formatter->endNode();
             return;
         }
+
+        auto &originalVertCRS = CRS::getPrivate()->originalVertCRS_;
+        if (!formatter->useESRIDialect() && originalVertCRS) {
+            formatter->startNode(io::WKTConstants::COMPD_CS, false);
+            formatter->addQuotedString(l_name + " + " +
+                                       originalVertCRS->nameStr());
+            projCRS2D->_exportToWKT(formatter);
+            originalVertCRS->_exportToWKT(formatter);
+            formatter->endNode();
+            return;
+        }
+
         io::FormattingException::Throw(
             "Projected 3D CRS can only be exported since WKT2:2019");
     }
@@ -4184,6 +4254,14 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties,
         auto comp1 = components[1].get();
         auto comp0Geog = dynamic_cast<const GeographicCRS *>(comp0);
         auto comp0Proj = dynamic_cast<const ProjectedCRS *>(comp0);
+        if (comp0Geog == nullptr && comp0Proj == nullptr) {
+            auto comp0Bound = dynamic_cast<const BoundCRS *>(comp0);
+            if (comp0Bound) {
+                const auto *baseCRS = comp0Bound->baseCRS().get();
+                comp0Geog = dynamic_cast<const GeographicCRS *>(baseCRS);
+                comp0Proj = dynamic_cast<const ProjectedCRS *>(baseCRS);
+            }
+        }
         auto comp1Geog = dynamic_cast<const GeographicCRS *>(comp1);
         if ((comp0Geog != nullptr || comp0Proj != nullptr) &&
             comp1Geog != nullptr) {
@@ -4201,6 +4279,22 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties,
                 "The 'vertical' geographic CRS is not equivalent to the "
                 "geographic CRS of the horizontal part");
         }
+
+        // Detect a COMPD_CS whose VERT_CS is for ellipoidal heights
+        auto comp1Vert =
+            util::nn_dynamic_pointer_cast<VerticalCRS>(components[1]);
+        if (comp1Vert != nullptr && comp1Vert->datum() &&
+            comp1Vert->datum()->getWKT1DatumType() == "2002") {
+            const auto &axis = comp1Vert->coordinateSystem()->axisList()[0];
+            if (axis->unit()._isEquivalentTo(
+                    common::UnitOfMeasure::METRE,
+                    util::IComparable::Criterion::EQUIVALENT) &&
+                &(axis->direction()) == &(cs::AxisDirection::UP)) {
+                return components[0]
+                    ->promoteTo3D(std::string(), dbContext)
+                    ->attachOriginalVertCRS(NN_NO_CHECK(comp1Vert));
+            }
+        }
     }
 
     return create(properties, components);
diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp
index 03abb0d7..33eb81fa 100644
--- a/src/iso19111/datum.cpp
+++ b/src/iso19111/datum.cpp
@@ -1788,6 +1788,9 @@ operator=(const RealizationMethod &other) {
 //! @cond Doxygen_Suppress
 struct VerticalReferenceFrame::Private {
     util::optional<RealizationMethod> realizationMethod_{};
+
+    // 2005 = CS_VD_GeoidModelDerived from OGC 01-009
+    std::string wkt1DatumType_{"2005"};
 };
 //! @endcond
 
@@ -1837,12 +1840,21 @@ VerticalReferenceFrameNNPtr VerticalReferenceFrame::create(
         realizationMethodIn));
     rf->setAnchor(anchor);
     rf->setProperties(properties);
+    properties.getStringValue("VERT_DATUM_TYPE", rf->d->wkt1DatumType_);
     return rf;
 }
 
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+const std::string &VerticalReferenceFrame::getWKT1DatumType() const {
+    return d->wkt1DatumType_;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 void VerticalReferenceFrame::_exportToWKT(
     io::WKTFormatter *formatter) const // throw(FormattingException)
 {
@@ -1861,7 +1873,7 @@ void VerticalReferenceFrame::_exportToWKT(
     if (isWKT2) {
         Datum::getPrivate()->exportAnchorDefinition(formatter);
     } else if (!formatter->useESRIDialect()) {
-        formatter->add(2005); // CS_VD_GeoidModelDerived from OGC 01-009
+        formatter->add(d->wkt1DatumType_);
         const auto &extension = formatter->getVDatumExtension();
         if (!extension.empty()) {
             formatter->startNode(io::WKTConstants::EXTENSION, false);
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 9d9fa4a3..ebec053a 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3944,9 +3944,18 @@ VerticalReferenceFrameNNPtr WKTParser::Private::buildVerticalReferenceFrame(
             modelName);
     }
 
-    // WKT1 VERT_DATUM has a datum type after the datum name that we ignore.
-    return VerticalReferenceFrame::create(buildProperties(node),
-                                          getAnchor(node));
+    // WKT1 VERT_DATUM has a datum type after the datum name
+    const auto *nodeP = node->GP();
+    const std::string &name(nodeP->value());
+    auto &props = buildProperties(node);
+    if (ci_equal(name, WKTConstants::VERT_DATUM)) {
+        const auto &children = nodeP->children();
+        if (children.size() >= 2) {
+            props.set("VERT_DATUM_TYPE", children[1]->GP()->value());
+        }
+    }
+
+    return VerticalReferenceFrame::create(props, getAnchor(node));
 }
 
 // ---------------------------------------------------------------------------