Merge pull request #2229 from rouault/fix_2228

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

6 files changed, 270 insertions, 20 deletions

diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 61844f79..7aa74c41 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -140,6 +140,9 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage,
 
     PROJ_INTERNAL CRSNNPtr allowNonConformantWKT1Export() const;
 
+    PROJ_INTERNAL CRSNNPtr
+    attachOriginalVertCRS(const VerticalCRSNNPtr &vertCRS) const;
+
     //! @endcond
 
   protected:
diff --git a/include/proj/datum.hpp b/include/proj/datum.hpp
index 6a0db1dc..3fc5f5d8 100644
--- a/include/proj/datum.hpp
+++ b/include/proj/datum.hpp
@@ -590,6 +590,8 @@ class PROJ_GCC_DLL VerticalReferenceFrame : public Datum {
     PROJ_INTERNAL void _exportToJSON(io::JSONFormatter *formatter)
         const override; // throw(FormattingException)
 
+    PROJ_INTERNAL const std::string &getWKT1DatumType() const;
+
     //! @endcond
 
   protected:
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));
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 7ad9cf4a..d62e5c59 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -2484,6 +2484,136 @@ TEST(wkt_parse, COMPD_CS_non_conformant_horizontal_plus_horizontal_as_in_LAS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse,
+     COMPD_CS_non_conformant_horizontal_TOWGS84_plus_horizontal_as_in_LAS) {
+
+    const auto wkt = "COMPD_CS[\"WGS 84 + WGS 84\",\n"
+                     "    GEOGCS[\"WGS 84\",\n"
+                     "        DATUM[\"WGS_1984\",\n"
+                     "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                     "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                     "            TOWGS84[0,0,0,0,0,0,0],\n"
+                     "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                     "        PRIMEM[\"Greenwich\",0,\n"
+                     "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                     "        UNIT[\"degree\",0.0174532925199433,\n"
+                     "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                     "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
+                     "    GEOGCS[\"WGS 84\",\n"
+                     "        DATUM[\"WGS_1984\",\n"
+                     "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                     "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                     "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                     "        PRIMEM[\"Greenwich\",0,\n"
+                     "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                     "        UNIT[\"degree\",0.0174532925199433,\n"
+                     "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                     "        AUTHORITY[\"EPSG\",\"4326\"]]]";
+    auto dbContext = DatabaseContext::create();
+    auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto baseCRS = nn_dynamic_pointer_cast<GeographicCRS>(crs->baseCRS());
+    ASSERT_TRUE(baseCRS != nullptr);
+    EXPECT_EQ(baseCRS->nameStr(), "WGS 84");
+    EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U);
+
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+                .get()),
+        wkt);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse,
+     COMPD_CS_horizontal_bound_geog_plus_vertical_ellipsoidal_height) {
+    // See https://github.com/OSGeo/PROJ/issues/2228
+    const char *wkt =
+        "COMPD_CS[\"NAD83 + Ellipsoid (Meters)\",\n"
+        "    GEOGCS[\"NAD83\",\n"
+        "        DATUM[\"North_American_Datum_1983\",\n"
+        "            SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                AUTHORITY[\"EPSG\",\"7019\"]],\n"
+        "            TOWGS84[0,0,0,0,0,0,0],\n"
+        "            AUTHORITY[\"EPSG\",\"6269\"]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+        "        UNIT[\"degree\",0.0174532925199433,\n"
+        "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+        "        AUTHORITY[\"EPSG\",\"4269\"]],\n"
+        "    VERT_CS[\"Ellipsoid (Meters)\",\n"
+        "        VERT_DATUM[\"Ellipsoid\",2002],\n"
+        "        UNIT[\"metre\",1,\n"
+        "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+        "        AXIS[\"Up\",UP]]]";
+    auto dbContext = DatabaseContext::create();
+    auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto baseCRS = nn_dynamic_pointer_cast<GeographicCRS>(crs->baseCRS());
+    ASSERT_TRUE(baseCRS != nullptr);
+    EXPECT_EQ(baseCRS->nameStr(), "NAD83");
+    EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U);
+
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+                .get()),
+        wkt);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse,
+     COMPD_CS_horizontal_projected_plus_vertical_ellipsoidal_height) {
+    // Variant of above
+    const char *wkt =
+        "COMPD_CS[\"WGS 84 / UTM zone 31N + Ellipsoid (Meters)\",\n"
+        "    PROJCS[\"WGS 84 / UTM zone 31N\",\n"
+        "        GEOGCS[\"WGS 84\",\n"
+        "            DATUM[\"WGS_1984\",\n"
+        "                SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+        "                    AUTHORITY[\"EPSG\",\"7030\"]],\n"
+        "                AUTHORITY[\"EPSG\",\"6326\"]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                AUTHORITY[\"EPSG\",\"8901\"]],\n"
+        "            UNIT[\"degree\",0.0174532925199433,\n"
+        "                AUTHORITY[\"EPSG\",\"9122\"]],\n"
+        "            AUTHORITY[\"EPSG\",\"4326\"]],\n"
+        "        PROJECTION[\"Transverse_Mercator\"],\n"
+        "        PARAMETER[\"latitude_of_origin\",0],\n"
+        "        PARAMETER[\"central_meridian\",3],\n"
+        "        PARAMETER[\"scale_factor\",0.9996],\n"
+        "        PARAMETER[\"false_easting\",500000],\n"
+        "        PARAMETER[\"false_northing\",0],\n"
+        "        UNIT[\"metre\",1,\n"
+        "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+        "        AXIS[\"Easting\",EAST],\n"
+        "        AXIS[\"Northing\",NORTH],\n"
+        "        AUTHORITY[\"EPSG\",\"32631\"]],\n"
+        "    VERT_CS[\"Ellipsoid (Meters)\",\n"
+        "        VERT_DATUM[\"Ellipsoid\",2002],\n"
+        "        UNIT[\"metre\",1,\n"
+        "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+        "        AXIS[\"Up\",UP]]]";
+    auto dbContext = DatabaseContext::create();
+    auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->nameStr(), "WGS 84 / UTM zone 31N");
+    EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U);
+
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+                .get()),
+        wkt);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, COORDINATEOPERATION) {
 
     std::string src_wkt;