Geographic 3D CRS: allow to export to WKT1:ESRI if only the GEOGCS is known (and thus extrapolating a VERTCS) (fixes #2757)

2 files changed, 63 insertions, 4 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 0d157da1..575f6e2b 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1850,11 +1850,39 @@ static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
     auto vertCRSList = authFactory->createVerticalCRSFromDatum(
         "ESRI", "from_geogdatum_" + *gdatum_ids[0]->codeSpace() + '_' +
                     gdatum_ids[0]->code());
-    if (vertCRSList.size() != 1) {
-        return false;
-    }
     self->demoteTo2D(std::string(), dbContext)->_exportToWKT(formatter);
-    vertCRSList.front()->_exportToWKT(formatter);
+    if (vertCRSList.size() == 1) {
+        vertCRSList.front()->_exportToWKT(formatter);
+    } else {
+        // This will not be recognized properly by ESRI software
+        // See https://github.com/OSGeo/PROJ/issues/2757
+
+        const auto &axisList = geodCRS->coordinateSystem()->axisList();
+        assert(axisList.size() == 3U);
+
+        formatter->startNode(io::WKTConstants::VERTCS, false);
+        auto vertcs_name = l_esri_name;
+        if (starts_with(vertcs_name.c_str(), "GCS_"))
+            vertcs_name = vertcs_name.substr(4);
+        formatter->addQuotedString(vertcs_name);
+
+        gdatum->_exportToWKT(formatter);
+
+        // Seems to be a constant value...
+        formatter->startNode(io::WKTConstants::PARAMETER, false);
+        formatter->addQuotedString("Vertical_Shift");
+        formatter->add(0.0);
+        formatter->endNode();
+
+        formatter->startNode(io::WKTConstants::PARAMETER, false);
+        formatter->addQuotedString("Direction");
+        formatter->add(
+            axisList[2]->direction() == cs::AxisDirection::UP ? 1.0 : -1.0);
+        formatter->endNode();
+
+        axisList[2]->unit()._exportToWKT(formatter);
+        formatter->endNode();
+    }
     return true;
 }
 
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index d55a8c6e..500ff4a4 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -787,6 +787,37 @@ TEST(crs,
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, implicit_compound_ESRI_104971_to_3D_as_WKT1_ESRI_with_database) {
+    auto dbContext = DatabaseContext::create();
+    auto factory = AuthorityFactory::create(dbContext, "ESRI");
+    auto crs = factory->createGeographicCRS("104971")->promoteTo3D(
+        std::string(), dbContext);
+    WKTFormatterNNPtr f(WKTFormatter::create(
+        WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
+    // Situation where there is no ESRI vertical CRS, but the GEOGCS does exist
+    // This will be only partly recognized by ESRI software.
+    // See https://github.com/OSGeo/PROJ/issues/2757
+    const char *wkt = "GEOGCS[\"Mars_2000_(Sphere)\","
+                      "DATUM[\"Mars_2000_(Sphere)\","
+                      "SPHEROID[\"Mars_2000_(Sphere)\",3396190.0,0.0]],"
+                      "PRIMEM[\"Reference_Meridian\",0.0],"
+                      "UNIT[\"Degree\",0.0174532925199433]],"
+                      "VERTCS[\"Mars_2000_(Sphere)\","
+                      "DATUM[\"Mars_2000_(Sphere)\","
+                      "SPHEROID[\"Mars_2000_(Sphere)\",3396190.0,0.0]],"
+                      "PARAMETER[\"Vertical_Shift\",0.0],"
+                      "PARAMETER[\"Direction\",1.0],"
+                      "UNIT[\"Meter\",1.0]]";
+    EXPECT_EQ(crs->exportToWKT(f.get()), wkt);
+
+    auto obj2 = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto crs2 = nn_dynamic_pointer_cast<GeographicCRS>(obj2);
+    ASSERT_TRUE(crs2 != nullptr);
+    EXPECT_EQ(crs2->coordinateSystem()->axisList().size(), 3U);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, IAU_1000_as_WKT2) {
     auto dbContext = DatabaseContext::create();
     auto factory = AuthorityFactory::create(dbContext, "IAU_2015");