PROJ4 string import: take into correctly non-metre unit when the string looks like the one for WGS 84 / Pseudo Mercator (fixes https://github.com/OSGeo/gdal/issues/2433) (#2174)

3 files changed, 109 insertions, 23 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index f61215d0..883e2f32 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -5821,6 +5821,24 @@ static bool createPROJ4WebMercator(const Conversion *conv,
         return false;
     }
 
+    std::string units("m");
+    auto targetCRS = conv->targetCRS();
+    auto targetProjCRS =
+        dynamic_cast<const crs::ProjectedCRS *>(targetCRS.get());
+    if (targetProjCRS) {
+        const auto &axisList = targetProjCRS->coordinateSystem()->axisList();
+        const auto &unit = axisList[0]->unit();
+        if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE,
+                                  util::IComparable::Criterion::EQUIVALENT)) {
+            auto projUnit = unit.exportToPROJString();
+            if (!projUnit.empty()) {
+                units = projUnit;
+            } else {
+                return false;
+            }
+        }
+    }
+
     formatter->addStep("merc");
     const double a = geogCRS->ellipsoid()->semiMajorAxis().getSIValue();
     formatter->addParam("a", a);
@@ -5830,7 +5848,7 @@ static bool createPROJ4WebMercator(const Conversion *conv,
     formatter->addParam("x_0", falseEasting);
     formatter->addParam("y_0", falseNorthing);
     formatter->addParam("k", 1.0);
-    formatter->addParam("units", "m");
+    formatter->addParam("units", units);
     formatter->addParam("nadgrids", "@null");
     formatter->addParam("wktext");
     formatter->addParam("no_defs");
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index aa676c3d..6ec0617c 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3650,13 +3650,13 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
 
 // ---------------------------------------------------------------------------
 
-static ProjectedCRSNNPtr createPseudoMercator(const PropertyMap &props) {
+static ProjectedCRSNNPtr createPseudoMercator(const PropertyMap &props,
+                                              const cs::CartesianCSNNPtr &cs) {
     auto conversion = Conversion::createPopularVisualisationPseudoMercator(
         PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), Angle(0),
         Angle(0), Length(0), Length(0));
-    return ProjectedCRS::create(
-        props, GeographicCRS::EPSG_4326, conversion,
-        CartesianCS::createEastingNorthing(UnitOfMeasure::METRE));
+    return ProjectedCRS::create(props, GeographicCRS::EPSG_4326, conversion,
+                                cs);
 }
 
 // ---------------------------------------------------------------------------
@@ -3682,6 +3682,15 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
 
     auto props = buildProperties(node);
 
+    auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
+    const auto &nodeValue = nodeP->value();
+    if (isNull(csNode) && !ci_equal(nodeValue, WKTConstants::PROJCS) &&
+        !ci_equal(nodeValue, WKTConstants::BASEPROJCRS)) {
+        ThrowMissing(WKTConstants::CS_);
+    }
+    auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
+    auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
+
     const std::string projCRSName = stripQuotes(nodeP->children()[0]);
     if (esriStyle_ && dbContext_) {
         // It is likely that the ESRI definition of EPSG:32661 (UPS North) &
@@ -3707,21 +3716,22 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
         }
     }
 
-    if (isNull(conversionNode) && hasWebMercPROJ4String(node, projectionNode)) {
+    if (isNull(conversionNode) && hasWebMercPROJ4String(node, projectionNode) &&
+        cartesianCS) {
         toWGS84Parameters_.clear();
-        return createPseudoMercator(props);
+        return createPseudoMercator(props, NN_NO_CHECK(cartesianCS));
     }
 
     // WGS_84_Pseudo_Mercator: Particular case for corrupted ESRI WKT generated
     // by older GDAL versions
     // https://trac.osgeo.org/gdal/changeset/30732
     // WGS_1984_Web_Mercator: deprecated ESRI:102113
-    if (metadata::Identifier::isEquivalentName(projCRSName.c_str(),
-                                               "WGS_84_Pseudo_Mercator") ||
-        metadata::Identifier::isEquivalentName(projCRSName.c_str(),
-                                               "WGS_1984_Web_Mercator")) {
+    if (cartesianCS && (metadata::Identifier::isEquivalentName(
+                            projCRSName.c_str(), "WGS_84_Pseudo_Mercator") ||
+                        metadata::Identifier::isEquivalentName(
+                            projCRSName.c_str(), "WGS_1984_Web_Mercator"))) {
         toWGS84Parameters_.clear();
-        return createPseudoMercator(props);
+        return createPseudoMercator(props, NN_NO_CHECK(cartesianCS));
     }
 
     auto linearUnit = buildUnitInSubNode(node, UnitOfMeasure::Type::LINEAR);
@@ -3733,15 +3743,6 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
             : buildProjection(baseGeodCRS, node, projectionNode, linearUnit,
                               angularUnit);
 
-    auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
-    const auto &nodeValue = nodeP->value();
-    if (isNull(csNode) && !ci_equal(nodeValue, WKTConstants::PROJCS) &&
-        !ci_equal(nodeValue, WKTConstants::BASEPROJCRS)) {
-        ThrowMissing(WKTConstants::CS_);
-    }
-    auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
-    auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
-
     // No explicit AXIS node ? (WKT1)
     if (isNull(nodeP->lookForChild(WKTConstants::AXIS))) {
         props.set("IMPLICIT_CS", true);
@@ -8597,6 +8598,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
 
     auto axisType = AxisType::REGULAR;
     bool bWebMercator = false;
+    std::string webMercatorName("WGS 84 / Pseudo-Mercator");
 
     if (step.name == "tmerc" &&
         ((getParamValue(step, "axis") == "wsu" && iAxisSwap < 0) ||
@@ -8685,6 +8687,11 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                 getAngularValue(getParamValue(step, "x_0")) == 0.0 &&
                 getAngularValue(getParamValue(step, "y_0")) == 0.0) {
                 bWebMercator = true;
+                if (hasParamValue(step, "units") &&
+                    getParamValue(step, "units") != "m") {
+                    webMercatorName +=
+                        " (unit " + getParamValue(step, "units") + ')';
+                }
             }
         } else if (hasParamValue(step, "lat_ts")) {
             mapping = getMapping(EPSG_CODE_METHOD_MERCATOR_VARIANT_B);
@@ -8968,8 +8975,8 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
 
     CRSNNPtr crs =
         bWebMercator
-            ? createPseudoMercator(props.set(IdentifiedObject::NAME_KEY,
-                                             "WGS 84 / Pseudo-Mercator"))
+            ? createPseudoMercator(
+                  props.set(IdentifiedObject::NAME_KEY, webMercatorName), cs)
             : ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs);
 
     if (!hasParamValue(step, "geoidgrids") &&
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index a2c75e04..8db48a83 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -8294,6 +8294,67 @@ TEST(io, projparse_merc_google_mercator) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_merc_google_mercator_non_metre_unit) {
+    auto projString =
+        "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 "
+        "+k=1 +units=ft +nadgrids=@null +no_defs +type=crs";
+    auto obj = PROJStringParser().createFromPROJString(projString);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->nameStr(), "WGS 84 / Pseudo-Mercator (unit ft)");
+    EXPECT_EQ(crs->derivingConversion()->method()->nameStr(),
+              "Popular Visualisation Pseudo Mercator");
+    EXPECT_EQ(
+        replaceAll(crs->exportToPROJString(PROJStringFormatter::create().get()),
+                   " +wktext", ""),
+        projString);
+
+    auto wkt = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get());
+    auto expected_wkt =
+        "PROJCRS[\"WGS 84 / Pseudo-Mercator (unit ft)\",\n"
+        "    BASEGEOGCRS[\"WGS 84\",\n"
+        "        DATUM[\"World Geodetic System 1984\",\n"
+        "            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                LENGTHUNIT[\"metre\",1]]],\n"
+        "        PRIMEM[\"Greenwich\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "        ID[\"EPSG\",4326]],\n"
+        "    CONVERSION[\"unnamed\",\n"
+        "        METHOD[\"Popular Visualisation Pseudo Mercator\",\n"
+        "            ID[\"EPSG\",1024]],\n"
+        "        PARAMETER[\"Latitude of natural origin\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "            ID[\"EPSG\",8801]],\n"
+        "        PARAMETER[\"Longitude of natural origin\",0,\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "            ID[\"EPSG\",8802]],\n"
+        "        PARAMETER[\"False easting\",0,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8806]],\n"
+        "        PARAMETER[\"False northing\",0,\n"
+        "            LENGTHUNIT[\"metre\",1],\n"
+        "            ID[\"EPSG\",8807]]],\n"
+        "    CS[Cartesian,2],\n"
+        "        AXIS[\"(E)\",east,\n"
+        "            ORDER[1],\n"
+        "            LENGTHUNIT[\"foot\",0.3048,\n"
+        "                ID[\"EPSG\",9002]]],\n"
+        "        AXIS[\"(N)\",north,\n"
+        "            ORDER[2],\n"
+        "            LENGTHUNIT[\"foot\",0.3048,\n"
+        "                ID[\"EPSG\",9002]]]]";
+    EXPECT_EQ(wkt, expected_wkt);
+
+    auto objFromWkt = WKTParser().createFromWKT(wkt);
+    auto crsFromWkt = nn_dynamic_pointer_cast<ProjectedCRS>(objFromWkt);
+    ASSERT_TRUE(crsFromWkt != nullptr);
+    EXPECT_TRUE(crs->isEquivalentTo(crsFromWkt.get(),
+                                    IComparable::Criterion::EQUIVALENT));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_merc_not_quite_google_mercator) {
     auto projString =
         "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=10 +x_0=0 +y_0=0 "