normalizeForVisualization(): address EPSG:32661 "WGS 84 / UPS North (N,E)" and EPSG:32761 "WGS 84 / UPS South (N,E)"

Fixes https://github.com/qgis/QGIS/issues/33077

2 files changed, 81 insertions, 7 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 49cc050f..ad4728b8 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -595,12 +595,44 @@ CRSNNPtr CRS::alterId(const std::string &authName,
 
 //! @cond Doxygen_Suppress
 
-static bool isAxisListNorthEast(
+static bool mustAxisOrderBeSwitchedForVisualizationInternal(
     const std::vector<cs::CoordinateSystemAxisNNPtr> &axisList) {
     const auto &dir0 = axisList[0]->direction();
     const auto &dir1 = axisList[1]->direction();
-    return (&dir0 == &cs::AxisDirection::NORTH &&
-            &dir1 == &cs::AxisDirection::EAST);
+    if (&dir0 == &cs::AxisDirection::NORTH &&
+        &dir1 == &cs::AxisDirection::EAST) {
+        return true;
+    }
+
+    // Address EPSG:32661 "WGS 84 / UPS North (N,E)"
+    if (&dir0 == &cs::AxisDirection::SOUTH &&
+        &dir1 == &cs::AxisDirection::SOUTH) {
+        const auto &meridian0 = axisList[0]->meridian();
+        const auto &meridian1 = axisList[1]->meridian();
+        return meridian0 != nullptr && meridian1 != nullptr &&
+               fabs(meridian0->longitude().convertToUnit(
+                        common::UnitOfMeasure::DEGREE) -
+                    180.0) < 1e-10 &&
+               fabs(meridian1->longitude().convertToUnit(
+                        common::UnitOfMeasure::DEGREE) -
+                    90.0) < 1e-10;
+    }
+
+    // Address EPSG:32761 "WGS 84 / UPS South (N,E)"
+    if (&dir0 == &cs::AxisDirection::NORTH &&
+        &dir1 == &cs::AxisDirection::NORTH) {
+        const auto &meridian0 = axisList[0]->meridian();
+        const auto &meridian1 = axisList[1]->meridian();
+        return meridian0 != nullptr && meridian1 != nullptr &&
+               fabs(meridian0->longitude().convertToUnit(
+                        common::UnitOfMeasure::DEGREE) -
+                    0.0) < 1e-10 &&
+               fabs(meridian1->longitude().convertToUnit(
+                        common::UnitOfMeasure::DEGREE) -
+                    90.0) < 1e-10;
+    }
+
+    return false;
 }
 // ---------------------------------------------------------------------------
 
@@ -616,12 +648,14 @@ bool CRS::mustAxisOrderBeSwitchedForVisualization() const {
 
     const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this);
     if (geogCRS) {
-        return isAxisListNorthEast(geogCRS->coordinateSystem()->axisList());
+        return mustAxisOrderBeSwitchedForVisualizationInternal(
+            geogCRS->coordinateSystem()->axisList());
     }
 
     const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this);
     if (projCRS) {
-        return isAxisListNorthEast(projCRS->coordinateSystem()->axisList());
+        return mustAxisOrderBeSwitchedForVisualizationInternal(
+            projCRS->coordinateSystem()->axisList());
     }
 
     return false;
@@ -655,7 +689,7 @@ CRSNNPtr CRS::normalizeForVisualization() const {
     const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this);
     if (geogCRS) {
         const auto &axisList = geogCRS->coordinateSystem()->axisList();
-        if (isAxisListNorthEast(axisList)) {
+        if (mustAxisOrderBeSwitchedForVisualizationInternal(axisList)) {
             auto cs = axisList.size() == 2
                           ? cs::EllipsoidalCS::create(util::PropertyMap(),
                                                       axisList[1], axisList[0])
@@ -670,7 +704,7 @@ CRSNNPtr CRS::normalizeForVisualization() const {
     const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this);
     if (projCRS) {
         const auto &axisList = projCRS->coordinateSystem()->axisList();
-        if (isAxisListNorthEast(axisList)) {
+        if (mustAxisOrderBeSwitchedForVisualizationInternal(axisList)) {
             auto cs =
                 axisList.size() == 2
                     ? cs::CartesianCS::create(util::PropertyMap(), axisList[1],
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index e896853f..33d67e0a 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -943,6 +943,16 @@ TEST(crs, EPSG_32661_projected_north_pole_north_east) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -964,6 +974,16 @@ TEST(crs, EPSG_5041_projected_north_pole_east_north) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -985,6 +1005,16 @@ TEST(crs, EPSG_32761_projected_south_pole_north_east) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------
@@ -1006,6 +1036,16 @@ TEST(crs, EPSG_5042_projected_south_pole_east_north) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               proj_string);
+
+    auto opNormalized = op->normalizeForVisualization();
+    auto proj_string_normalized =
+        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step "
+        "+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 "
+        "+ellps=WGS84";
+    EXPECT_EQ(
+        opNormalized->exportToPROJString(PROJStringFormatter::create().get()),
+        proj_string_normalized);
 }
 
 // ---------------------------------------------------------------------------