import from WKT1: fix axis direction when importing from Polar Stereographic variant B

2 files changed, 86 insertions, 35 deletions

diff --git a/src/io.cpp b/src/io.cpp
index b5da2b8c..f25ece48 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -3460,6 +3460,37 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
     }
 
     if (isNull(csNode) && node->countChildrenOfName(WKTConstants::AXIS) == 0) {
+
+        const auto createEastingSouthNorthingSouth = [&linearUnit]() {
+            return CartesianCS::create(
+                PropertyMap(),
+                CoordinateSystemAxis::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            AxisName::Easting),
+                    AxisAbbreviation::E, AxisDirection::SOUTH, linearUnit,
+                    Meridian::create(Angle(90, UnitOfMeasure::DEGREE))),
+                CoordinateSystemAxis::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            AxisName::Northing),
+                    AxisAbbreviation::N, AxisDirection::SOUTH, linearUnit,
+                    Meridian::create(Angle(180, UnitOfMeasure::DEGREE))));
+        };
+
+        const auto createEastingNorthNorthingNorth = [&linearUnit]() {
+            return CartesianCS::create(
+                PropertyMap(),
+                CoordinateSystemAxis::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            AxisName::Easting),
+                    AxisAbbreviation::E, AxisDirection::NORTH, linearUnit,
+                    Meridian::create(Angle(90, UnitOfMeasure::DEGREE))),
+                CoordinateSystemAxis::create(
+                    util::PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                            AxisName::Northing),
+                    AxisAbbreviation::N, AxisDirection::NORTH, linearUnit,
+                    Meridian::create(Angle(0, UnitOfMeasure::DEGREE))));
+        };
+
         const auto methodCode = conversion->method()->getEPSGCode();
         // Krovak south oriented ?
         if (methodCode == EPSG_CODE_METHOD_KROVAK) {
@@ -3475,49 +3506,29 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
                                                 AxisName::Westing),
                         emptyString, AxisDirection::WEST, linearUnit))
                     .as_nullable();
-        } else if (esriStyle_ &&
-                   methodCode ==
-                       EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_A) {
+        } else if (methodCode ==
+                   EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_A) {
             // 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.
+            // Same for WKT1_GDAL
             const double lat0 = conversion->parameterValueNumeric(
                 EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN,
                 common::UnitOfMeasure::DEGREE);
             if (std::fabs(lat0 - 90) < 1e-10) {
-                cartesianCS =
-                    CartesianCS::create(
-                        PropertyMap(),
-                        CoordinateSystemAxis::create(
-                            util::PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                                    AxisName::Easting),
-                            AxisAbbreviation::E, AxisDirection::SOUTH,
-                            linearUnit,
-                            Meridian::create(Angle(90, UnitOfMeasure::DEGREE))),
-                        CoordinateSystemAxis::create(
-                            util::PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                                    AxisName::Northing),
-                            AxisAbbreviation::N, AxisDirection::SOUTH,
-                            linearUnit, Meridian::create(
-                                            Angle(180, UnitOfMeasure::DEGREE))))
-                        .as_nullable();
+                cartesianCS = createEastingSouthNorthingSouth().as_nullable();
             } else if (std::fabs(lat0 - -90) < 1e-10) {
-                cartesianCS =
-                    CartesianCS::create(
-                        PropertyMap(),
-                        CoordinateSystemAxis::create(
-                            util::PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                                    AxisName::Easting),
-                            AxisAbbreviation::E, AxisDirection::NORTH,
-                            linearUnit,
-                            Meridian::create(Angle(90, UnitOfMeasure::DEGREE))),
-                        CoordinateSystemAxis::create(
-                            util::PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                                    AxisName::Northing),
-                            AxisAbbreviation::N, AxisDirection::NORTH,
-                            linearUnit,
-                            Meridian::create(Angle(0, UnitOfMeasure::DEGREE))))
-                        .as_nullable();
+                cartesianCS = createEastingNorthNorthingNorth().as_nullable();
+            }
+        } else if (methodCode ==
+                   EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_B) {
+            const double lat_ts = conversion->parameterValueNumeric(
+                EPSG_CODE_PARAMETER_LATITUDE_STD_PARALLEL,
+                common::UnitOfMeasure::DEGREE);
+            if (lat_ts > 0) {
+                cartesianCS = createEastingSouthNorthingSouth().as_nullable();
+            } else if (lat_ts < 0) {
+                cartesianCS = createEastingNorthNorthingNorth().as_nullable();
             }
         }
     }
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 1f2ec323..32d72823 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -1356,6 +1356,46 @@ TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_70) {
     auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=2 "
                               "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs";
     EXPECT_EQ(projString, expectedPROJString);
+
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[0]->nameStr(), "Easting");
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[0]->direction(),
+              AxisDirection::SOUTH);
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[1]->nameStr(), "Northing");
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[1]->direction(),
+              AxisDirection::SOUTH);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_minus_70) {
+    auto wkt = "PROJCS[\"unknown\",\n"
+               "    GEOGCS[\"unknown\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]]],\n"
+               "    PROJECTION[\"Polar_Stereographic\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",-70],\n"
+               "    PARAMETER[\"central_meridian\",2],\n"
+               "    PARAMETER[\"false_easting\",3],\n"
+               "    PARAMETER[\"false_northing\",4],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[0]->nameStr(), "Easting");
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[0]->direction(),
+              AxisDirection::NORTH);
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[1]->nameStr(), "Northing");
+    EXPECT_EQ(crs->coordinateSystem()->axisList()[1]->direction(),
+              AxisDirection::NORTH);
 }
 
 // ---------------------------------------------------------------------------