Merge pull request #2902 from rouault/issue_2757

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

3 files changed, 218 insertions, 15 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 6d213bc6..575f6e2b 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1828,15 +1828,16 @@ static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
         return false;
     }
     const auto l_datum = geodCRS->datumNonNull(formatter->databaseContext());
-    auto l_alias = dbContext->getAliasFromOfficialName(
+    auto l_esri_name = dbContext->getAliasFromOfficialName(
         l_datum->nameStr(), "geodetic_datum", "ESRI");
-    if (l_alias.empty()) {
-        return false;
+    if (l_esri_name.empty()) {
+        l_esri_name = l_datum->nameStr();
     }
     auto authFactory =
         io::AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string());
     auto list = authFactory->createObjectsFromName(
-        l_alias, {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
+        l_esri_name,
+        {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
         false /* approximate=false*/);
     if (list.empty()) {
         return false;
@@ -1849,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;
 }
 
@@ -1983,6 +2012,18 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
                     aliasFound = true;
                 }
             }
+            if (!aliasFound && dbContext) {
+                auto authFactory = io::AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext), "ESRI");
+                aliasFound =
+                    authFactory
+                        ->createObjectsFromName(
+                            l_name,
+                            {io::AuthorityFactory::ObjectType::GEODETIC_CRS},
+                            false // approximateMatch
+                            )
+                        .size() == 1;
+            }
             if (!aliasFound) {
                 l_name = io::WKTFormatter::morphNameToESRI(l_name);
                 if (!starts_with(l_name, "GCS_")) {
@@ -3266,6 +3307,18 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
                 aliasFound = true;
             }
         }
+        if (!aliasFound && dbContext) {
+            auto authFactory =
+                io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "ESRI");
+            aliasFound =
+                authFactory
+                    ->createObjectsFromName(
+                        l_name,
+                        {io::AuthorityFactory::ObjectType::VERTICAL_CRS},
+                        false // approximateMatch
+                        )
+                    .size() == 1;
+        }
         if (!aliasFound) {
             l_name = io::WKTFormatter::morphNameToESRI(l_name);
         }
@@ -3966,10 +4019,24 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             "Projected 3D CRS can only be exported since WKT2:2019");
     }
 
-    std::string l_alias;
+    std::string l_esri_name;
     if (formatter->useESRIDialect() && dbContext) {
-        l_alias = dbContext->getAliasFromOfficialName(l_name, "projected_crs",
-                                                      "ESRI");
+        l_esri_name = dbContext->getAliasFromOfficialName(
+            l_name, "projected_crs", "ESRI");
+        if (l_esri_name.empty()) {
+            auto authFactory =
+                io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "ESRI");
+            const bool found =
+                authFactory
+                    ->createObjectsFromName(
+                        l_name,
+                        {io::AuthorityFactory::ObjectType::PROJECTED_CRS},
+                        false // approximateMatch
+                        )
+                    .size() == 1;
+            if (found)
+                l_esri_name = l_name;
+        }
     }
 
     if (!isWKT2 && formatter->useESRIDialect() && !l_identifiers.empty() &&
@@ -3991,12 +4058,12 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             }
         } catch (const std::exception &) {
         }
-    } else if (!isWKT2 && formatter->useESRIDialect() && !l_alias.empty()) {
+    } else if (!isWKT2 && formatter->useESRIDialect() && !l_esri_name.empty()) {
         try {
             auto res =
                 io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "ESRI")
                     ->createObjectsFromName(
-                        l_alias,
+                        l_esri_name,
                         {io::AuthorityFactory::ObjectType::PROJECTED_CRS},
                         false);
             if (res.size() == 1) {
@@ -4075,10 +4142,10 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
                          !l_identifiers.empty());
 
     if (formatter->useESRIDialect()) {
-        if (l_alias.empty()) {
+        if (l_esri_name.empty()) {
             l_name = io::WKTFormatter::morphNameToESRI(l_name);
         } else {
-            l_name = l_alias;
+            l_name = l_esri_name;
         }
     }
     if (!isWKT2 && !formatter->useESRIDialect() && isDeprecated()) {
diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp
index ebef94a2..7c76061e 100644
--- a/src/iso19111/datum.cpp
+++ b/src/iso19111/datum.cpp
@@ -363,6 +363,18 @@ void PrimeMeridian::_exportToWKT(
                     aliasFound = true;
                 }
             }
+            if (!aliasFound && dbContext) {
+                auto authFactory = io::AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext), "ESRI");
+                aliasFound =
+                    authFactory
+                        ->createObjectsFromName(
+                            l_name,
+                            {io::AuthorityFactory::ObjectType::PRIME_MERIDIAN},
+                            false // approximateMatch
+                            )
+                        .size() == 1;
+            }
             if (!aliasFound) {
                 l_name = io::WKTFormatter::morphNameToESRI(l_name);
             }
@@ -820,6 +832,18 @@ void Ellipsoid::_exportToWKT(
                             aliasFound = true;
                         }
                     }
+                    if (!aliasFound && dbContext) {
+                        auto authFactory = io::AuthorityFactory::create(
+                            NN_NO_CHECK(dbContext), "ESRI");
+                        aliasFound = authFactory
+                                         ->createObjectsFromName(
+                                             l_name,
+                                             {io::AuthorityFactory::ObjectType::
+                                                  ELLIPSOID},
+                                             false // approximateMatch
+                                             )
+                                         .size() == 1;
+                    }
                     if (!aliasFound) {
                         l_name = io::WKTFormatter::morphNameToESRI(l_name);
                     }
@@ -1247,6 +1271,18 @@ void GeodeticReferenceFrame::_exportToWKT(
                         }
                     }
                 }
+                if (!aliasFound && dbContext) {
+                    auto authFactory = io::AuthorityFactory::create(
+                        NN_NO_CHECK(dbContext), "ESRI");
+                    aliasFound = authFactory
+                                     ->createObjectsFromName(
+                                         l_name,
+                                         {io::AuthorityFactory::ObjectType::
+                                              GEODETIC_REFERENCE_FRAME},
+                                         false // approximateMatch
+                                         )
+                                     .size() == 1;
+                }
                 if (!aliasFound) {
                     l_name = io::WKTFormatter::morphNameToESRI(l_name);
                     if (!starts_with(l_name, "D_")) {
@@ -1985,6 +2021,18 @@ void VerticalReferenceFrame::_exportToWKT(
                     aliasFound = true;
                 }
             }
+            if (!aliasFound && dbContext) {
+                auto authFactory = io::AuthorityFactory::create(
+                    NN_NO_CHECK(dbContext), "ESRI");
+                aliasFound = authFactory
+                                 ->createObjectsFromName(
+                                     l_name,
+                                     {io::AuthorityFactory::ObjectType::
+                                          VERTICAL_REFERENCE_FRAME},
+                                     false // approximateMatch
+                                     )
+                                 .size() == 1;
+            }
             if (!aliasFound) {
                 l_name = io::WKTFormatter::morphNameToESRI(l_name);
             }
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 987ffe1b..500ff4a4 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -747,6 +747,77 @@ TEST(crs, EPSG_4268_geogcrs_deprecated_as_WKT1_GDAL) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, ESRI_104971_as_WKT1_ESRI_with_database) {
+    auto dbContext = DatabaseContext::create();
+    auto factory = AuthorityFactory::create(dbContext, "ESRI");
+    auto crs = factory->createCoordinateReferenceSystem("104971");
+    WKTFormatterNNPtr f(WKTFormatter::create(
+        WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
+    // Check that the _(Sphere) suffix is preserved
+    EXPECT_EQ(crs->exportToWKT(f.get()),
+              "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]]");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(crs,
+     implicit_compound_ESRI_104024_plus_115844_as_WKT1_ESRI_with_database) {
+    auto dbContext = DatabaseContext::create();
+    auto obj = createFromUserInput("ESRI:104024+115844", dbContext);
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U);
+    WKTFormatterNNPtr f(WKTFormatter::create(
+        WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
+    // Situation where there is no EPSG official name
+    EXPECT_EQ(crs->exportToWKT(f.get()),
+              "GEOGCS[\"California_SRS_Epoch_2017.50_(NAD83)\","
+              "DATUM[\"California_SRS_Epoch_2017.50_(NAD83)\","
+              "SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],"
+              "PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],"
+              "VERTCS[\"California_SRS_Epoch_2017.50_(NAD83)\","
+              "DATUM[\"California_SRS_Epoch_2017.50_(NAD83)\","
+              "SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],"
+              "PARAMETER[\"Vertical_Shift\",0.0],"
+              "PARAMETER[\"Direction\",1.0],UNIT[\"Meter\",1.0]]");
+}
+
+// ---------------------------------------------------------------------------
+
+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");
@@ -3877,6 +3948,23 @@ TEST(crs, verticalCRS_down_as_WKT1_ESRI) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, verticalCRS_ESRI_115834_as_WKT1_ESRI_with_database) {
+    auto dbContext = DatabaseContext::create();
+    auto factory = AuthorityFactory::create(dbContext, "ESRI");
+    auto crs = factory->createCoordinateReferenceSystem("115834");
+    WKTFormatterNNPtr f(WKTFormatter::create(
+        WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
+    // Check that the parentheses in the VERTCS and DATUM names are preserved
+    EXPECT_EQ(crs->exportToWKT(f.get()),
+              "VERTCS[\"NAD83(CSRS)v5\","
+              "DATUM[\"North_American_Datum_of_1983_(CSRS)_version_5\","
+              "SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],"
+              "PARAMETER[\"Vertical_Shift\",0.0],"
+              "PARAMETER[\"Direction\",1.0],UNIT[\"Meter\",1.0]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, verticalCRS_identify_db) {
     auto dbContext = DatabaseContext::create();
     auto factory = AuthorityFactory::create(dbContext, "EPSG");