proj_normalize_for_visualization(): fix when there are coordinate operation alternatives

2 files changed, 117 insertions, 0 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 6a2a1ae0..85421fa0 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -30,6 +30,7 @@
 #define FROM_PROJ_CPP
 #endif
 
+#include <algorithm>
 #include <cassert>
 #include <cstdarg>
 #include <cstring>
@@ -151,6 +152,7 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) {
     }
     auto pj = pj_new();
     if (pj) {
+        pj->ctx = ctx;
         pj->descr = "ISO-19111 object";
         pj->iso_obj = objIn;
     }
@@ -6793,6 +6795,56 @@ int proj_cs_get_axis_info(PJ_CONTEXT *ctx, const PJ *cs, int index,
  * nullptr in case of error
  */
 PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) {
+
+    if (!obj->alternativeCoordinateOperations.empty()) {
+        try {
+            auto pjNew = pj_new();
+            pjNew->ctx = ctx;
+            for (const auto &alt : obj->alternativeCoordinateOperations) {
+                auto co = dynamic_cast<const CoordinateOperation *>(
+                    alt.pj->iso_obj.get());
+                if (co) {
+                    double minxSrc = alt.minxSrc;
+                    double minySrc = alt.minySrc;
+                    double maxxSrc = alt.maxxSrc;
+                    double maxySrc = alt.maxySrc;
+                    double minxDst = alt.minxDst;
+                    double minyDst = alt.minyDst;
+                    double maxxDst = alt.maxxDst;
+                    double maxyDst = alt.maxyDst;
+
+                    auto l_sourceCRS = co->sourceCRS();
+                    auto l_targetCRS = co->targetCRS();
+                    if (l_sourceCRS && l_targetCRS) {
+                        const bool swapSource =
+                            l_sourceCRS
+                                ->mustAxisOrderBeSwitchedForVisualization();
+                        if (swapSource) {
+                            std::swap(minxSrc, minySrc);
+                            std::swap(maxxSrc, maxySrc);
+                        }
+                        const bool swapTarget =
+                            l_targetCRS
+                                ->mustAxisOrderBeSwitchedForVisualization();
+                        if (swapTarget) {
+                            std::swap(minxDst, minyDst);
+                            std::swap(maxxDst, maxyDst);
+                        }
+                    }
+                    pjNew->alternativeCoordinateOperations.emplace_back(
+                        minxSrc, minySrc, maxxSrc, maxySrc, minxDst, minyDst,
+                        maxxDst, maxyDst,
+                        pj_obj_create(ctx, co->normalizeForVisualization()),
+                        co->nameStr());
+                }
+            }
+            return pjNew;
+        } catch (const std::exception &e) {
+            proj_log_debug(ctx, __FUNCTION__, e.what());
+            return nullptr;
+        }
+    }
+
     auto co = dynamic_cast<const CoordinateOperation *>(obj->iso_obj.get());
     if (!co) {
         proj_log_error(ctx, __FUNCTION__, "Object is not a CoordinateOperation "
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 12e98d7e..bea3eaa1 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -3247,4 +3247,69 @@ TEST_F(CApi, proj_normalize_for_visualization) {
               "+step +proj=utm +zone=31 +ellps=WGS84");
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_normalize_for_visualization_with_alternatives) {
+
+    auto P = proj_create_crs_to_crs(m_ctxt, "EPSG:4326", "EPSG:3003", nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper keeper_Pnormalized(Pnormalized);
+    ASSERT_NE(Pnormalized, nullptr);
+
+    {
+        PJ_COORD c;
+        // Approximately Roma
+        c.lpz.lam = 12.5;
+        c.lpz.phi = 42;
+        c.lpz.z = 0;
+        c = proj_trans(Pnormalized, PJ_FWD, c);
+        EXPECT_NEAR(c.xy.x, 1789912.46264783037, 1e-8);
+        EXPECT_NEAR(c.xy.y, 4655716.25402576849, 1e-8);
+        auto projstr = proj_pj_info(Pnormalized).definition;
+        ASSERT_NE(projstr, nullptr);
+        EXPECT_EQ(std::string(projstr),
+                  "proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad "
+                  "step proj=push v_3 step proj=cart ellps=WGS84 "
+                  "step inv proj=helmert x=-104.1 y=-49.1 z=-9.9 rx=0.971 "
+                  "ry=-2.917 rz=0.714 s=-11.68 convention=position_vector "
+                  "step inv proj=cart ellps=intl step proj=pop v_3 "
+                  "step proj=tmerc lat_0=0 lon_0=9 k=0.9996 x_0=1500000 "
+                  "y_0=0 ellps=intl");
+    }
+
+    {
+        PJ_COORD c;
+        // Approximately Roma
+        c.xyz.x = 1789912.46264783037;
+        c.xyz.y = 4655716.25402576849;
+        c.xyz.z = 0;
+        c = proj_trans(Pnormalized, PJ_INV, c);
+        EXPECT_NEAR(c.lp.lam, 12.5, 1e-8);
+        EXPECT_NEAR(c.lp.phi, 42, 1e-8);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_normalize_for_visualization_with_alternatives_reverse) {
+
+    auto P = proj_create_crs_to_crs(m_ctxt, "EPSG:3003", "EPSG:4326", nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto Pnormalized = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper keeper_Pnormalized(Pnormalized);
+    ASSERT_NE(Pnormalized, nullptr);
+
+    PJ_COORD c;
+    // Approximately Roma
+    c.xyz.x = 1789912.46264783037;
+    c.xyz.y = 4655716.25402576849;
+    c.xyz.z = 0;
+    c = proj_trans(Pnormalized, PJ_FWD, c);
+    EXPECT_NEAR(c.lp.lam, 12.5, 1e-8);
+    EXPECT_NEAR(c.lp.phi, 42, 1e-8);
+}
+
 } // namespace