PRIMEM WKT handling: fixes on import for 'sexagesimal DMS' or from WKT1:GDAL/ESRI when GEOGCS UNIT != Degree; morph to ESRI the PRIMEM name on export

4 files changed, 159 insertions, 4 deletions

diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp
index 60244779..8e1b8257 100644
--- a/src/iso19111/datum.cpp
+++ b/src/iso19111/datum.cpp
@@ -351,6 +351,23 @@ void PrimeMeridian::_exportToWKT(
     if (!(isWKT2 && formatter->primeMeridianOmittedIfGreenwich() &&
           l_name == "Greenwich")) {
         formatter->startNode(io::WKTConstants::PRIMEM, !identifiers().empty());
+
+        if (formatter->useESRIDialect()) {
+            bool aliasFound = false;
+            const auto &dbContext = formatter->databaseContext();
+            if (dbContext) {
+                auto l_alias = dbContext->getAliasFromOfficialName(
+                    l_name, "prime_meridian", "ESRI");
+                if (!l_alias.empty()) {
+                    l_name = l_alias;
+                    aliasFound = true;
+                }
+            }
+            if (!aliasFound) {
+                l_name = io::WKTFormatter::morphNameToESRI(l_name);
+            }
+        }
+
         formatter->addQuotedString(l_name);
         const auto &l_long = longitude();
         if (formatter->primeMeridianInDegree()) {
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 867a0d53..52d46d27 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -1978,14 +1978,80 @@ PrimeMeridianNNPtr WKTParser::Private::buildPrimeMeridian(
     try {
         double angleValue = asDouble(children[1]);
 
-        // Correct for GDAL WKT1 departure
+        // Correct for GDAL WKT1 and WKT1-ESRI departure
         if (name == "Paris" && std::fabs(angleValue - 2.33722917) < 1e-8 &&
-            unit == UnitOfMeasure::GRAD) {
+            unit._isEquivalentTo(UnitOfMeasure::GRAD,
+                                 util::IComparable::Criterion::EQUIVALENT)) {
             angleValue = 2.5969213;
+        } else {
+            static const struct {
+                const char *name;
+                int deg;
+                int min;
+                double sec;
+            } primeMeridiansDMS[] = {
+                {"Lisbon", -9, 7, 54.862},  {"Bogota", -74, 4, 51.3},
+                {"Madrid", -3, 41, 14.55},  {"Rome", 12, 27, 8.4},
+                {"Bern", 7, 26, 22.5},      {"Jakarta", 106, 48, 27.79},
+                {"Ferro", -17, 40, 0},      {"Brussels", 4, 22, 4.71},
+                {"Stockholm", 18, 3, 29.8}, {"Athens", 23, 42, 58.815},
+                {"Oslo", 10, 43, 22.5},     {"Paris RGS", 2, 20, 13.95},
+                {"Paris_RGS", 2, 20, 13.95}};
+
+            // Current epsg.org output may use the EPSG:9110 "sexagesimal DMS"
+            // unit and a DD.MMSSsss value, but this will likely be changed to
+            // use decimal degree.
+            // Or WKT1 may for example use the Paris RGS decimal degree value
+            // but with a GEOGCS with UNIT["Grad"]
+            for (const auto &pmDef : primeMeridiansDMS) {
+                if (name == pmDef.name) {
+                    double dmsAsDecimalValue =
+                        (pmDef.deg >= 0 ? 1 : -1) *
+                        (std::abs(pmDef.deg) + pmDef.min / 100. +
+                         pmDef.sec / 10000.);
+                    double dmsAsDecimalDegreeValue =
+                        (pmDef.deg >= 0 ? 1 : -1) *
+                        (std::abs(pmDef.deg) + pmDef.min / 60. +
+                         pmDef.sec / 3600.);
+                    if (std::fabs(angleValue - dmsAsDecimalValue) < 1e-8 ||
+                        std::fabs(angleValue - dmsAsDecimalDegreeValue) <
+                            1e-8) {
+                        angleValue = dmsAsDecimalDegreeValue;
+                        unit = UnitOfMeasure::DEGREE;
+                    }
+                    break;
+                }
+            }
+        }
+
+        auto &properties = buildProperties(node);
+        if (dbContext_ && esriStyle_) {
+            std::string outTableName;
+            std::string codeFromAlias;
+            std::string authNameFromAlias;
+            auto authFactory = AuthorityFactory::create(NN_NO_CHECK(dbContext_),
+                                                        std::string());
+            auto officialName = authFactory->getOfficialNameFromAlias(
+                name, "prime_meridian", "ESRI", false, outTableName,
+                authNameFromAlias, codeFromAlias);
+            if (!officialName.empty()) {
+                properties.set(IdentifiedObject::NAME_KEY, officialName);
+                if (!authNameFromAlias.empty()) {
+                    auto identifiers = ArrayOfBaseObject::create();
+                    identifiers->add(Identifier::create(
+                        codeFromAlias,
+                        PropertyMap()
+                            .set(Identifier::CODESPACE_KEY, authNameFromAlias)
+                            .set(Identifier::AUTHORITY_KEY,
+                                 authNameFromAlias)));
+                    properties.set(IdentifiedObject::IDENTIFIERS_KEY,
+                                   identifiers);
+                }
+            }
         }
 
         Angle angle(angleValue, unit);
-        return PrimeMeridian::create(buildProperties(node), angle);
+        return PrimeMeridian::create(properties, angle);
     } catch (const std::exception &e) {
         throw buildRethrow(__FUNCTION__, e);
     }
@@ -2737,6 +2803,9 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
         throw ParsingException("Missing DATUM or ENSEMBLE node");
     }
 
+    // Do that now so that esriStyle_ can be set before buildPrimeMeridian()
+    auto props = buildProperties(node);
+
     auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC);
 
     auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
@@ -2774,7 +2843,6 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
         angularUnit = primeMeridian->longitude().unit();
     }
 
-    auto props = buildProperties(node);
     addExtensionProj4ToProp(nodeP, props);
 
     // No explicit AXIS node ? (WKT1)
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index e470a621..af172b4c 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -414,6 +414,20 @@ TEST(crs, EPSG_4326_as_WKT1_ESRI_with_database) {
 
 // ---------------------------------------------------------------------------
 
+TEST(crs, EPSG_4901_as_WKT1_ESRI_with_PRIMEM_unit_name_morphing) {
+    auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto crs = factory->createCoordinateReferenceSystem("4901");
+    WKTFormatterNNPtr f(WKTFormatter::create(
+        WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
+    EXPECT_EQ(crs->exportToWKT(f.get()),
+              "GEOGCS[\"GCS_ATF_Paris\",DATUM[\"D_ATF\","
+              "SPHEROID[\"Plessis_1817\",6376523.0,308.64]],"
+              "PRIMEM[\"Paris_RGS\",2.33720833333333],"
+              "UNIT[\"Grad\",0.0157079632679489]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, EPSG_4326_as_WKT1_ESRI_without_database) {
     auto crs = GeographicCRS::EPSG_4326;
     WKTFormatterNNPtr f(
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 81894fb0..f0e221d8 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -488,6 +488,62 @@ TEST(wkt_parse, wkt1_EPSG_4807_grad_mess) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_esri_EPSG_4901_grad) {
+    auto obj =
+        WKTParser()
+            .attachDatabaseContext(DatabaseContext::create())
+            .createFromWKT("GEOGCS[\"GCS_ATF_Paris\",DATUM[\"D_ATF\","
+                           "SPHEROID[\"Plessis_1817\",6376523.0,308.64]],"
+                           "PRIMEM[\"Paris_RGS\",2.33720833333333],"
+                           "UNIT[\"Grad\",0.0157079632679489]]");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto datum = crs->datum();
+    auto primem = datum->primeMeridian();
+    EXPECT_EQ(primem->nameStr(), "Paris RGS");
+    // The PRIMEM is really in degree
+    EXPECT_EQ(primem->longitude().unit(), UnitOfMeasure::DEGREE);
+    EXPECT_NEAR(primem->longitude().value(), 2.33720833333333, 1e-14);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt2_epsg_org_EPSG_4901_PRIMEM_weird_sexagesimal_DMS) {
+    // Current epsg.org output may use the EPSG:9110 "sexagesimal DMS"
+    // unit and a DD.MMSSsss value, but this will likely be changed to
+    // use decimal degree.
+    auto obj = WKTParser().createFromWKT(
+        "GEOGCRS[\"ATF (Paris)\","
+        "  DATUM[\"Ancienne Triangulation Francaise (Paris)\","
+        "    ELLIPSOID[\"Plessis 1817\",6376523,308.64,"
+        "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],"
+        "      ID[\"EPSG\",7027]],"
+        "    ID[\"EPSG\",6901]],"
+        "  PRIMEM[\"Paris RGS\",2.201395,"
+        "    ANGLEUNIT[\"sexagesimal DMS\",1,ID[\"EPSG\",9110]],"
+        "    ID[\"EPSG\",8914]],"
+        "  CS[ellipsoidal,2,"
+        "    ID[\"EPSG\",6403]],"
+        "  AXIS[\"Geodetic latitude (Lat)\",north,"
+        "    ORDER[1]],"
+        "  AXIS[\"Geodetic longitude (Lon)\",east,"
+        "    ORDER[2]],"
+        "  ANGLEUNIT[\"grad\",0.015707963267949,ID[\"EPSG\",9105]],"
+        "  USAGE[SCOPE[\"Geodesy.\"],AREA[\"France - mainland onshore.\"],"
+        "  BBOX[42.33,-4.87,51.14,8.23]],"
+        "ID[\"EPSG\",4901]]");
+    auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    auto datum = crs->datum();
+    auto primem = datum->primeMeridian();
+    EXPECT_EQ(primem->longitude().unit(), UnitOfMeasure::DEGREE);
+    EXPECT_NEAR(primem->longitude().value(), 2.33720833333333, 1e-14);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt1_geographic_old_datum_name_from_EPSG_code) {
     auto wkt =
         "GEOGCS[\"S-JTSK (Ferro)\",\n"