Support Projected 3D with method name suffixed with ' (3D)', and be more lax on projection parameter names

3 files changed, 139 insertions, 6 deletions

diff --git a/src/coordinateoperation.cpp b/src/coordinateoperation.cpp
index 98a2d830..442a9b78 100644
--- a/src/coordinateoperation.cpp
+++ b/src/coordinateoperation.cpp
@@ -159,6 +159,10 @@ static std::set<std::string> buildSetEquivalentParameters() {
          EPSG_NAME_PARAMETER_NORTHING_FALSE_ORIGIN,
          EPSG_NAME_PARAMETER_NORTHING_PROJECTION_CENTRE, nullptr},
 
+        {EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN, WKT1_SCALE_FACTOR,
+         EPSG_NAME_PARAMETER_SCALE_FACTOR_INITIAL_LINE,
+         EPSG_NAME_PARAMETER_SCALE_FACTOR_PSEUDO_STANDARD_PARALLEL, nullptr},
+
         {WKT1_LATITUDE_OF_ORIGIN, WKT1_LATITUDE_OF_CENTER,
          EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN,
          EPSG_NAME_PARAMETER_LATITUDE_FALSE_ORIGIN,
@@ -826,7 +830,10 @@ OperationMethodNNPtr OperationMethod::create(
 int OperationMethod::getEPSGCode() PROJ_CONST_DEFN {
     int epsg_code = IdentifiedObject::getEPSGCode();
     if (epsg_code == 0) {
-        const auto &l_name = nameStr();
+        auto l_name = nameStr();
+        if (ends_with(l_name, " (3D)")) {
+            l_name.resize(l_name.size() - strlen(" (3D)"));
+        }
         for (const auto &tuple : methodNameCodes) {
             if (metadata::Identifier::isEquivalentName(l_name.c_str(),
                                                        tuple.name)) {
@@ -1268,6 +1275,14 @@ int OperationParameter::getEPSGCode() PROJ_CONST_DEFN {
                 return tuple.epsg_code;
             }
         }
+        if (metadata::Identifier::isEquivalentName(l_name.c_str(),
+                                                   "Latitude of origin")) {
+            return EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN;
+        }
+        if (metadata::Identifier::isEquivalentName(l_name.c_str(),
+                                                   "Scale factor")) {
+            return EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN;
+        }
     }
     return epsg_code;
 }
@@ -1350,13 +1365,35 @@ static const ParameterValuePtr nullParameterValue;
 const ParameterValuePtr &
 SingleOperation::parameterValue(const std::string &paramName,
                                 int epsg_code) const noexcept {
+    if (epsg_code) {
+        for (const auto &genOpParamvalue : parameterValues()) {
+            auto opParamvalue = dynamic_cast<const OperationParameterValue *>(
+                genOpParamvalue.get());
+            if (opParamvalue) {
+                const auto &parameter = opParamvalue->parameter();
+                if (parameter->getEPSGCode() == epsg_code) {
+                    return opParamvalue->parameterValue();
+                }
+            }
+        }
+    }
+    for (const auto &genOpParamvalue : parameterValues()) {
+        auto opParamvalue = dynamic_cast<const OperationParameterValue *>(
+            genOpParamvalue.get());
+        if (opParamvalue) {
+            const auto &parameter = opParamvalue->parameter();
+            if (metadata::Identifier::isEquivalentName(
+                    paramName.c_str(), parameter->nameStr().c_str())) {
+                return opParamvalue->parameterValue();
+            }
+        }
+    }
     for (const auto &genOpParamvalue : parameterValues()) {
         auto opParamvalue = dynamic_cast<const OperationParameterValue *>(
             genOpParamvalue.get());
         if (opParamvalue) {
             const auto &parameter = opParamvalue->parameter();
-            if ((epsg_code != 0 && parameter->getEPSGCode() == epsg_code) ||
-                ci_equal(paramName, parameter->nameStr())) {
+            if (areEquivalentParameters(paramName, parameter->nameStr())) {
                 return opParamvalue->parameterValue();
             }
         }
@@ -5739,8 +5776,10 @@ bool Conversion::isUTM(int &zone, bool &north) const {
                                   UTM_LATITUDE_OF_NATURAL_ORIGIN) < 1e-10) {
                         bLatitudeNatOriginUTM = true;
                     } else if (
-                        epsg_code ==
-                            EPSG_CODE_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN &&
+                        (epsg_code ==
+                             EPSG_CODE_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN ||
+                         epsg_code ==
+                             EPSG_CODE_PARAMETER_LONGITUDE_OF_ORIGIN) &&
                         measure.unit()._isEquivalentTo(
                             common::UnitOfMeasure::DEGREE,
                             util::IComparable::Criterion::EQUIVALENT)) {
diff --git a/src/crs.cpp b/src/crs.cpp
index d873e81b..55539e18 100644
--- a/src/crs.cpp
+++ b/src/crs.cpp
@@ -1060,7 +1060,7 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     if (formatter->useESRIDialect()) {
         if (axisList.size() != 2) {
             io::FormattingException::Throw(
-                "Only export of Geographic 2D CRS is supported in ESRI_WKT1");
+                "Only export of Geographic 2D CRS is supported in WKT1_ESRI");
         }
 
         if (l_name == "WGS 84") {
@@ -2590,6 +2590,10 @@ 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->use2018Keywords())) {
+        io::FormattingException::Throw(
+            "Projected 3D CRS can only be exported since WKT2:2018");
+    }
 
     const auto exportAxis = [&l_coordinateSystem, &axisList, &formatter]() {
         const auto oldAxisOutputRule = formatter->outputAxis();
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 1ed855f3..466ce80b 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -1617,6 +1617,96 @@ TEST(wkt_parse, wkt2_2018_simplified_projected) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt2_2018_projected_3D) {
+    auto wkt =
+        "PROJCRS[\"WGS 84 (G1762) / UTM zone 31N 3D\","
+        "    BASEGEOGCRS[\"WGS 84\","
+        "        DATUM[\"World Geodetic System of 1984 (G1762)\","
+        "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
+        "           LENGTHUNIT[\"metre\",1.0]]]],"
+        "    CONVERSION[\"Some conversion 3D\","
+        "        METHOD[\"Transverse Mercator (3D)\"],"
+        "        PARAMETER[\"Latitude of origin\",0.0,"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
+        "        PARAMETER[\"Longitude of origin\",3.0,"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
+        "        PARAMETER[\"Scale factor\",1,SCALEUNIT[\"unity\",1.0]],"
+        "        PARAMETER[\"False easting\",0.0,"
+        "           LENGTHUNIT[\"metre\",1.0]],"
+        "        PARAMETER[\"False northing\",0.0,LENGTHUNIT[\"metre\",1.0]]],"
+        "        CS[Cartesian,3],"
+        "            AXIS[\"(E)\",east,ORDER[1]],"
+        "            AXIS[\"(N)\",north,ORDER[2]],"
+        "            AXIS[\"ellipsoidal height (h)\",up,ORDER[3]],"
+        "            LENGTHUNIT[\"metre\",1.0]"
+        "]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=tmerc +lat_0=0 +lon_0=3 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 "
+        "+units=m +no_defs");
+
+    EXPECT_THROW(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT2_2015).get()),
+        FormattingException);
+
+    EXPECT_NO_THROW(crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT2_2018).get()));
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt2_2018_projected_utm_3D) {
+    // Example from WKT2:2018
+    auto wkt =
+        "PROJCRS[\"WGS 84 (G1762) / UTM zone 31N 3D\","
+        "    BASEGEOGCRS[\"WGS 84\","
+        "        DATUM[\"World Geodetic System of 1984 (G1762)\","
+        "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
+        "           LENGTHUNIT[\"metre\",1.0]]]],"
+        "    CONVERSION[\"UTM zone 31N 3D\","
+        "        METHOD[\"Transverse Mercator (3D)\"],"
+        "        PARAMETER[\"Latitude of origin\",0.0,"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
+        "        PARAMETER[\"Longitude of origin\",3.0,"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
+        "        PARAMETER[\"Scale factor\",0.9996,SCALEUNIT[\"unity\",1.0]],"
+        "        PARAMETER[\"False easting\",500000.0,"
+        "           LENGTHUNIT[\"metre\",1.0]],"
+        "        PARAMETER[\"False northing\",0.0,LENGTHUNIT[\"metre\",1.0]]],"
+        "        CS[Cartesian,3],"
+        "            AXIS[\"(E)\",east,ORDER[1]],"
+        "            AXIS[\"(N)\",north,ORDER[2]],"
+        "            AXIS[\"ellipsoidal height (h)\",up,ORDER[3]],"
+        "            LENGTHUNIT[\"metre\",1.0]"
+        "]";
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=utm +zone=31 +ellps=WGS84 +units=m +no_defs");
+
+    EXPECT_THROW(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT2_2015).get()),
+        FormattingException);
+
+    EXPECT_NO_THROW(crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT2_2018).get()));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, projected_angular_unit_from_primem) {
     auto obj = WKTParser().createFromWKT(
         "PROJCRS[\"NTF (Paris) / Lambert Nord France\",\n"