From 048cd320b4e63b9378cf6332c8bff11f03a6c39b Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 13 Dec 2018 15:51:12 +0100
Subject: import from WKT1 and PROJ string: tune CS for Polar LAEA

---
 src/io.cpp            | 15 ++++++++++++++-
 test/unit/test_io.cpp | 30 ++++++++++++++++++++++++++++++
 2 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/src/io.cpp b/src/io.cpp
index 1700b1e3..71d029ea 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -3477,7 +3477,9 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
                         emptyString, AxisDirection::WEST, linearUnit))
                     .as_nullable();
         } else if (methodCode ==
-                   EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_A) {
+                       EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_A ||
+                   methodCode ==
+                       EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA) {
             // It is likely that the ESRI definition of EPSG:32661 (UPS North) &
             // EPSG:32761 (UPS South) uses the easting-northing order, instead
             // of the EPSG northing-easting order.
@@ -6585,6 +6587,17 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
         } else {
             mapping = getMapping(PROJ_WKT2_NAME_METHOD_STEREOGRAPHIC);
         }
+    } else if (step.name == "laea") {
+        if (hasParamValue(step, "lat_0") &&
+            std::fabs(std::fabs(getAngularValue(getParamValue(step, "lat_0"))) -
+                      90.0) < 1e-10) {
+            const double lat_0 = getAngularValue(getParamValue(step, "lat_0"));
+            if (lat_0 > 0) {
+                axisType = AxisType::NORTH_POLE;
+            } else {
+                axisType = AxisType::SOUTH_POLE;
+            }
+        }
     }
 
     UnitOfMeasure unit = buildUnit(step, "units", "to_meter");
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 32d72823..bd230b45 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -6943,6 +6943,36 @@ TEST(io, projparse_utm_south) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_laea_north_pole) {
+    auto obj = PROJStringParser().createFromPROJString("+proj=laea +lat_0=90");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    WKTFormatterNNPtr f(WKTFormatter::create());
+    f->simulCurNodeHasId();
+    f->setMultiLine(false);
+    crs->exportToWKT(f.get());
+    auto wkt = f->toString();
+    EXPECT_TRUE(wkt.find("AXIS[\"(E)\",south") != std::string::npos) << wkt;
+    EXPECT_TRUE(wkt.find("AXIS[\"(N)\",south") != std::string::npos) << wkt;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_laea_south_pole) {
+    auto obj = PROJStringParser().createFromPROJString("+proj=laea +lat_0=-90");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    WKTFormatterNNPtr f(WKTFormatter::create());
+    f->simulCurNodeHasId();
+    f->setMultiLine(false);
+    crs->exportToWKT(f.get());
+    auto wkt = f->toString();
+    EXPECT_TRUE(wkt.find("AXIS[\"(E)\",north") != std::string::npos) << wkt;
+    EXPECT_TRUE(wkt.find("AXIS[\"(N)\",north") != std::string::npos) << wkt;
+}
+
+// ---------------------------------------------------------------------------
+
 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");
-- 
cgit v1.2.3