Merge pull request #2167 from rouault/backport_7_0_pr2166

[Backport 7.0] Ingestion of WKT1_GDAL: correctly map 'Cylindrical_Equal_Area'

5 files changed, 116 insertions, 32 deletions

diff --git a/include/proj/internal/esri_projection_mappings.hpp b/include/proj/internal/esri_projection_mappings.hpp
index 707f9c47..fdd3ef30 100644
--- a/include/proj/internal/esri_projection_mappings.hpp
+++ b/include/proj/internal/esri_projection_mappings.hpp
@@ -898,9 +898,8 @@ static const ESRIMethodMapping esriMappings[] = {
      EPSG_NAME_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA,
      EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA,
      paramsESRI_Lambert_Azimuthal_Equal_Area},
-    {"Cylindrical_Equal_Area",
-     EPSG_NAME_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL,
-     EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL,
+    {"Cylindrical_Equal_Area", EPSG_NAME_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA,
+     EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA,
      paramsESRI_Cylindrical_Equal_Area},
     {"Hotine_Oblique_Mercator_Two_Point_Center",
      PROJ_WKT2_NAME_METHOD_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN, 0,
diff --git a/scripts/build_esri_projection_mapping.py b/scripts/build_esri_projection_mapping.py
index 84eba426..dd183d54 100644
--- a/scripts/build_esri_projection_mapping.py
+++ b/scripts/build_esri_projection_mapping.py
@@ -330,7 +330,7 @@ config_str = """
         - Latitude_Of_Origin: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
 
 - Cylindrical_Equal_Area:
-    WKT2_name: EPSG_NAME_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL
+    WKT2_name: EPSG_NAME_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA
     Params:
         - False_Easting: EPSG_NAME_PARAMETER_FALSE_EASTING
         - False_Northing: EPSG_NAME_PARAMETER_FALSE_NORTHING
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d038b5f9..f8f824a9 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -1797,6 +1797,14 @@ bool SingleOperation::_isEquivalentTo(const util::IComparable *other,
                  EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA &&
              methodEPSGCode ==
                  EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA_SPHERICAL) ||
+            (methodEPSGCode ==
+                 EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA &&
+             otherMethodEPSGCode ==
+                 EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL) ||
+            (otherMethodEPSGCode ==
+                 EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA &&
+             methodEPSGCode ==
+                 EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL) ||
             (methodEPSGCode == EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL &&
              otherMethodEPSGCode ==
                  EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL_SPHERICAL) ||
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 9787bd30..de1c67c3 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -1298,19 +1298,22 @@ struct WKTParser::Private {
     static std::string projectionGetParameter(const WKTNodeNNPtr &projCRSNode,
                                               const char *paramName);
 
-    ConversionNNPtr buildProjection(const WKTNodeNNPtr &projCRSNode,
+    ConversionNNPtr buildProjection(const GeodeticCRSNNPtr &baseGeodCRS,
+                                    const WKTNodeNNPtr &projCRSNode,
                                     const WKTNodeNNPtr &projectionNode,
                                     const UnitOfMeasure &defaultLinearUnit,
                                     const UnitOfMeasure &defaultAngularUnit);
 
     ConversionNNPtr
-    buildProjectionStandard(const WKTNodeNNPtr &projCRSNode,
+    buildProjectionStandard(const GeodeticCRSNNPtr &baseGeodCRS,
+                            const WKTNodeNNPtr &projCRSNode,
                             const WKTNodeNNPtr &projectionNode,
                             const UnitOfMeasure &defaultLinearUnit,
                             const UnitOfMeasure &defaultAngularUnit);
 
     ConversionNNPtr
-    buildProjectionFromESRI(const WKTNodeNNPtr &projCRSNode,
+    buildProjectionFromESRI(const GeodeticCRSNNPtr &baseGeodCRS,
+                            const WKTNodeNNPtr &projCRSNode,
                             const WKTNodeNNPtr &projectionNode,
                             const UnitOfMeasure &defaultLinearUnit,
                             const UnitOfMeasure &defaultAngularUnit);
@@ -3138,9 +3141,40 @@ bool WKTParser::Private::hasWebMercPROJ4String(
 
 // ---------------------------------------------------------------------------
 
+static const MethodMapping *
+selectSphericalOrEllipsoidal(const MethodMapping *mapping,
+                             const GeodeticCRSNNPtr &baseGeodCRS) {
+    if (mapping->epsg_code ==
+            EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL ||
+        mapping->epsg_code == EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA) {
+        mapping = getMapping(
+            baseGeodCRS->ellipsoid()->isSphere()
+                ? EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL
+                : EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA);
+    } else if (mapping->epsg_code ==
+                   EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA_SPHERICAL ||
+               mapping->epsg_code ==
+                   EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA) {
+        mapping = getMapping(
+            baseGeodCRS->ellipsoid()->isSphere()
+                ? EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA_SPHERICAL
+                : EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA);
+    } else if (mapping->epsg_code ==
+                   EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL_SPHERICAL ||
+               mapping->epsg_code == EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL) {
+        mapping =
+            getMapping(baseGeodCRS->ellipsoid()->isSphere()
+                           ? EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL_SPHERICAL
+                           : EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL);
+    }
+    return mapping;
+}
+
+// ---------------------------------------------------------------------------
+
 ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
-    const WKTNodeNNPtr &projCRSNode, const WKTNodeNNPtr &projectionNode,
-    const UnitOfMeasure &defaultLinearUnit,
+    const GeodeticCRSNNPtr &baseGeodCRS, const WKTNodeNNPtr &projCRSNode,
+    const WKTNodeNNPtr &projectionNode, const UnitOfMeasure &defaultLinearUnit,
     const UnitOfMeasure &defaultAngularUnit) {
     const std::string esriProjectionName =
         stripQuotes(projectionNode->GP()->children()[0]);
@@ -3150,7 +3184,7 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
     // on the parameters / their values
     const auto esriMappings = getMappingsFromESRI(esriProjectionName);
     if (esriMappings.empty()) {
-        return buildProjectionStandard(projCRSNode, projectionNode,
+        return buildProjectionStandard(baseGeodCRS, projCRSNode, projectionNode,
                                        defaultLinearUnit, defaultAngularUnit);
     }
 
@@ -3230,6 +3264,8 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
     }
     assert(wkt2_mapping);
 
+    wkt2_mapping = selectSphericalOrEllipsoidal(wkt2_mapping, baseGeodCRS);
+
     PropertyMap propertiesMethod;
     propertiesMethod.set(IdentifiedObject::NAME_KEY, wkt2_mapping->wkt2_name);
     if (wkt2_mapping->epsg_code != 0) {
@@ -3317,19 +3353,18 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
 
 // ---------------------------------------------------------------------------
 
-ConversionNNPtr
-WKTParser::Private::buildProjection(const WKTNodeNNPtr &projCRSNode,
-                                    const WKTNodeNNPtr &projectionNode,
-                                    const UnitOfMeasure &defaultLinearUnit,
-                                    const UnitOfMeasure &defaultAngularUnit) {
+ConversionNNPtr WKTParser::Private::buildProjection(
+    const GeodeticCRSNNPtr &baseGeodCRS, const WKTNodeNNPtr &projCRSNode,
+    const WKTNodeNNPtr &projectionNode, const UnitOfMeasure &defaultLinearUnit,
+    const UnitOfMeasure &defaultAngularUnit) {
     if (projectionNode->GP()->childrenSize() == 0) {
         ThrowNotEnoughChildren(WKTConstants::PROJECTION);
     }
     if (esriStyle_) {
-        return buildProjectionFromESRI(projCRSNode, projectionNode,
+        return buildProjectionFromESRI(baseGeodCRS, projCRSNode, projectionNode,
                                        defaultLinearUnit, defaultAngularUnit);
     }
-    return buildProjectionStandard(projCRSNode, projectionNode,
+    return buildProjectionStandard(baseGeodCRS, projCRSNode, projectionNode,
                                    defaultLinearUnit, defaultAngularUnit);
 }
 
@@ -3354,8 +3389,8 @@ WKTParser::Private::projectionGetParameter(const WKTNodeNNPtr &projCRSNode,
 // ---------------------------------------------------------------------------
 
 ConversionNNPtr WKTParser::Private::buildProjectionStandard(
-    const WKTNodeNNPtr &projCRSNode, const WKTNodeNNPtr &projectionNode,
-    const UnitOfMeasure &defaultLinearUnit,
+    const GeodeticCRSNNPtr &baseGeodCRS, const WKTNodeNNPtr &projCRSNode,
+    const WKTNodeNNPtr &projectionNode, const UnitOfMeasure &defaultLinearUnit,
     const UnitOfMeasure &defaultAngularUnit) {
     std::string wkt1ProjectionName =
         stripQuotes(projectionNode->GP()->children()[0]);
@@ -3473,6 +3508,9 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard(
     std::string projectionName(wkt1ProjectionName);
     const MethodMapping *mapping =
         tryToIdentifyWKT1Method ? getMappingFromWKT1(projectionName) : nullptr;
+    if (mapping) {
+        mapping = selectSphericalOrEllipsoidal(mapping, baseGeodCRS);
+    }
 
     // For Krovak, we need to look at axis to decide between the Krovak and
     // Krovak East-North Oriented methods
@@ -3725,7 +3763,8 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
     auto conversion =
         !isNull(conversionNode)
             ? buildConversion(conversionNode, linearUnit, angularUnit)
-            : buildProjection(node, projectionNode, linearUnit, angularUnit);
+            : buildProjection(baseGeodCRS, node, projectionNode, linearUnit,
+                              angularUnit);
 
     // No explicit AXIS node ? (WKT1)
     if (isNull(nodeP->lookForChild(WKTConstants::AXIS))) {
@@ -8563,6 +8602,9 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     auto &step = steps_[iStep];
     auto mappings = getMappingsFromPROJName(step.name);
     const MethodMapping *mapping = mappings.empty() ? nullptr : mappings[0];
+    if (mapping) {
+        mapping = selectSphericalOrEllipsoidal(mapping, geogCRS);
+    }
 
     assert(isProjectedStep(step.name));
     assert(iUnitConvert < 0 ||
@@ -8608,8 +8650,6 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
         }
     } else if (step.name == "aeqd" && hasParamValue(step, "guam")) {
         mapping = getMapping(EPSG_CODE_METHOD_GUAM_PROJECTION);
-    } else if (step.name == "cea" && !geogCRS->ellipsoid()->isSphere()) {
-        mapping = getMapping(EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA);
     } else if (step.name == "geos" && getParamValue(step, "sweep") == "x") {
         mapping =
             getMapping(PROJ_WKT2_NAME_METHOD_GEOSTATIONARY_SATELLITE_SWEEP_X);
@@ -8715,15 +8755,6 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                 axisType = AxisType::SOUTH_POLE;
             }
         }
-        if (geogCRS->ellipsoid()->isSphere()) {
-            mapping = getMapping(
-                EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA_SPHERICAL);
-        }
-    } else if (step.name == "eqc") {
-        if (geogCRS->ellipsoid()->isSphere()) {
-            mapping =
-                getMapping(EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL_SPHERICAL);
-        }
     }
 
     UnitOfMeasure unit = buildUnit(step, "units", "to_meter");
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index ee828e62..fa4245ba 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -4658,7 +4658,7 @@ static const struct {
       {"False_Northing", 2},
       {"Central_Meridian", 3},
       {"Standard_Parallel_1", 4}},
-     "Lambert Cylindrical Equal Area (Spherical)",
+     "Lambert Cylindrical Equal Area",
      {
          {"Latitude of 1st standard parallel", 4},
          {"Longitude of natural origin", 3},
@@ -7969,6 +7969,30 @@ TEST(io, projparse_aeqd_guam) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_cea_spherical) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=cea +R=6371228 +type=crs");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->derivingConversion()->method()->getEPSGCode(),
+              EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL);
+
+    auto crs2 = ProjectedCRS::create(
+        PropertyMap(), crs->baseCRS(),
+        Conversion::createLambertCylindricalEqualArea(
+            PropertyMap(), Angle(0), Angle(0), Length(0), Length(0)),
+        crs->coordinateSystem());
+    EXPECT_EQ(crs2->derivingConversion()->method()->getEPSGCode(),
+              EPSG_CODE_METHOD_LAMBERT_CYLINDRICAL_EQUAL_AREA);
+
+    EXPECT_TRUE(
+        crs->isEquivalentTo(crs2.get(), IComparable::Criterion::EQUIVALENT));
+    EXPECT_TRUE(
+        crs2->isEquivalentTo(crs.get(), IComparable::Criterion::EQUIVALENT));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_cea_ellipsoidal) {
     auto obj = PROJStringParser().createFromPROJString(
         "+proj=cea +ellps=GRS80 +type=crs");
@@ -8492,6 +8516,17 @@ TEST(io, projparse_laea_spherical) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_laea_ellipsoidal) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=laea +ellps=WGS84 +type=crs");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->derivingConversion()->method()->getEPSGCode(),
+              EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_eqc_spherical) {
     auto obj = PROJStringParser().createFromPROJString(
         "+proj=eqc +R=6371228 +type=crs");
@@ -8516,6 +8551,17 @@ TEST(io, projparse_eqc_spherical) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_eqc_ellipsoidal) {
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=eqc +ellps=WGS84 +type=crs");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->derivingConversion()->method()->getEPSGCode(),
+              EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_non_earth_ellipsoid) {
     std::string input("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1 +units=m "
                       "+no_defs +type=crs");