Optimize pipelines of planetary CRS (geocentric latitude, west-positive longitude)

5 files changed, 155 insertions, 7 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 598cedad..a12c6aa0 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -2073,8 +2073,16 @@ void GeodeticCRS::_exportToPROJString(
         addGeocentricUnitConversionIntoPROJString(formatter);
     } else if (isSphericalPlanetocentric()) {
         if (!formatter->getCRSExport()) {
-            formatter->addStep("geoc");
-            addDatumInfoToPROJString(formatter);
+
+            if (!formatter->omitProjLongLatIfPossible() ||
+                primeMeridian()->longitude().getSIValue() != 0.0 ||
+                !ellipsoid()->isSphere() ||
+                !formatter->getTOWGS84Parameters().empty() ||
+                !formatter->getHDatumExtension().empty()) {
+                formatter->addStep("geoc");
+                addDatumInfoToPROJString(formatter);
+            }
+
             addAngularUnitConvertAndAxisSwap(formatter);
         } else {
             io::FormattingException::Throw(
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index d3928a86..b317213a 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -7509,13 +7509,22 @@ const std::string &PROJStringFormatter::toString() const {
 
         const auto paramCount = step.paramValues.size();
 
-        // axisswap order=2,1 is its own inverse
+        // axisswap order=2,1 (or 1,-2) is its own inverse
         if (step.name == "axisswap" && paramCount == 1 &&
-            step.paramValues[0].equals("order", "2,1")) {
+            (step.paramValues[0].equals("order", "2,1") ||
+             step.paramValues[0].equals("order", "1,-2"))) {
             step.inverted = false;
             continue;
         }
 
+        // axisswap inv order=2,-1 ==> axisswap order -2,1
+        if (step.name == "axisswap" && paramCount == 1 &&
+            step.paramValues[0].equals("order", "2,-1")) {
+            step.inverted = false;
+            step.paramValues[0] = Step::KeyValue("order", "-2,1");
+            continue;
+        }
+
         // handle unitconvert inverse
         if (step.name == "unitconvert" && paramCount == 2 &&
             step.paramValues[0].keyEquals("xy_in") &&
@@ -7790,6 +7799,25 @@ const std::string &PROJStringFormatter::toString() const {
                 continue;
             }
 
+            // axisswap order=2,1 followed by axisswap order=2,-1 is
+            // equivalent to axisswap order=1,-2
+            // Same for axisswap order=-2,1 followed by axisswap order=2,1
+            if (curStep.name == "axisswap" && prevStep.name == "axisswap" &&
+                curStepParamCount == 1 && prevStepParamCount == 1 &&
+                ((prevStep.paramValues[0].equals("order", "2,1") &&
+                  !curStep.inverted &&
+                  curStep.paramValues[0].equals("order", "2,-1")) ||
+                 (prevStep.paramValues[0].equals("order", "-2,1") &&
+                  !prevStep.inverted &&
+                  curStep.paramValues[0].equals("order", "2,1")))) {
+
+                prevStep.inverted = false;
+                prevStep.paramValues[0] = Step::KeyValue("order", "1,-2");
+                // Delete this iter
+                iterCur = steps.erase(iterCur);
+                continue;
+            }
+
             // axisswap order=2,1, unitconvert, axisswap order=2,1 -> can
             // suppress axisswap
             if (std::next(iterCur) != steps.end() &&
diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
index f7bb8ec8..3db62a7d 100644
--- a/src/iso19111/operation/conversion.cpp
+++ b/src/iso19111/operation/conversion.cpp
@@ -3534,11 +3534,13 @@ void Conversion::_exportToPROJString(
             if ((isSrcGeocentricLat && isTargetGeographic) ||
                 (isSrcGeographic && isTargetGeocentricLat)) {
 
+                formatter->setOmitProjLongLatIfPossible(true);
                 formatter->startInversion();
                 sourceCRSGeod->_exportToPROJString(formatter);
                 formatter->stopInversion();
 
                 targetCRSGeod->_exportToPROJString(formatter);
+                formatter->setOmitProjLongLatIfPossible(false);
 
                 return;
             }
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index ed059fe3..f9a54f3c 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -7852,6 +7852,56 @@ TEST(io, projstringformatter_axisswap_unitconvert_axisswap) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projstringformatter_axisswap_one_minus_two_inv) {
+    auto fmt = PROJStringFormatter::create();
+    fmt->ingestPROJString(
+        "+proj=pipeline +step +inv +proj=axisswap +order=1,-2");
+    EXPECT_EQ(fmt->toString(), "+proj=axisswap +order=1,-2");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projstringformatter_axisswap_two_one_followed_two_minus_one) {
+    auto fmt = PROJStringFormatter::create();
+    fmt->ingestPROJString("+proj=pipeline "
+                          "+step +proj=axisswap +order=2,1 "
+                          "+step +proj=axisswap +order=2,-1");
+    EXPECT_EQ(fmt->toString(), "+proj=axisswap +order=1,-2");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projstringformatter_axisswap_minus_two_one_followed_two_one) {
+    auto fmt = PROJStringFormatter::create();
+    fmt->ingestPROJString("+proj=pipeline "
+                          "+step +proj=axisswap +order=-2,1 "
+                          "+step +proj=axisswap +order=2,1");
+    EXPECT_EQ(fmt->toString(), "+proj=axisswap +order=1,-2");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projstringformatter_unmodified) {
+    const char *const strs[] = {"+proj=pipeline "
+                                "+step +proj=axisswap +order=2,-1 "
+                                "+step +proj=axisswap +order=2,1",
+
+                                "+proj=pipeline "
+                                "+step +proj=axisswap +order=2,1 "
+                                "+step +proj=axisswap +order=-2,1",
+
+                                "+proj=pipeline "
+                                "+step +inv +proj=axisswap +order=-2,1 "
+                                "+step +proj=axisswap +order=2,1"};
+    for (const char *str : strs) {
+        auto fmt = PROJStringFormatter::create();
+        fmt->ingestPROJString(str);
+        EXPECT_EQ(fmt->toString(), str);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projstringformatter_optim_hgridshift_vgridshift_hgridshift_inv) {
     // Nominal case
     {
@@ -8786,9 +8836,12 @@ TEST(io, projparse_longlat_axisswap) {
                     op->exportToPROJString(PROJStringFormatter::create().get()),
                     (atoi(order1) == 2 && atoi(order2) == 1)
                         ? "+proj=noop"
-                        : "+proj=pipeline +step +proj=axisswap +order=2,1 "
-                          "+step +proj=axisswap +order=" +
-                              std::string(order1) + "," + order2);
+                        : (atoi(order1) == 2 && atoi(order2) == -1)
+                              ? "+proj=axisswap +order=1,-2"
+                              : "+proj=pipeline +step +proj=axisswap "
+                                "+order=2,1 "
+                                "+step +proj=axisswap +order=" +
+                                    std::string(order1) + "," + order2);
             }
         }
     }
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index 2ee5a53c..31b23287 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -6953,3 +6953,60 @@ TEST(operation,
               "+step +inv +proj=geoc +ellps=WGS84 "
               "+step +proj=utm +zone=11 +ellps=WGS84");
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(operation,
+     createOperation_spherical_ocentric_spherical_to_ellipsoidal_north_west) {
+    auto objSrc = WKTParser().createFromWKT(
+        "GEODCRS[\"Mercury (2015) - Sphere / Ocentric\",\n"
+        "    DATUM[\"Mercury (2015) - Sphere\",\n"
+        "        ELLIPSOID[\"Mercury (2015) - Sphere\",2440530,0,\n"
+        "            LENGTHUNIT[\"metre\",1]]],\n"
+        "    PRIMEM[\"Reference Meridian\",0,\n"
+        "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "    CS[spherical,2],\n"
+        "        AXIS[\"planetocentric latitude (U)\",north,\n"
+        "            ORDER[1],\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "        AXIS[\"planetocentric longitude (V)\",east,\n"
+        "            ORDER[2],\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]]]");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest = WKTParser().createFromWKT(
+        "GEOGCRS[\"Mercury (2015) / Ographic\",\n"
+        "    DATUM[\"Mercury (2015)\",\n"
+        "        ELLIPSOID[\"Mercury (2015)\",2440530,1075.12334801762,\n"
+        "            LENGTHUNIT[\"metre\",1]]],\n"
+        "    PRIMEM[\"Reference Meridian\",0,\n"
+        "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "    CS[ellipsoidal,2],\n"
+        "        AXIS[\"latitude\",north,\n"
+        "            ORDER[1],\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+        "        AXIS[\"longitude\",west,\n"
+        "            ORDER[2],\n"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]]]"
+
+    );
+    auto dest = nn_dynamic_pointer_cast<CRS>(objDest);
+    ASSERT_TRUE(dest != nullptr);
+
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest));
+        ASSERT_TRUE(op != nullptr);
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=axisswap +order=1,-2");
+    }
+
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(dest), NN_CHECK_ASSERT(src));
+        ASSERT_TRUE(op != nullptr);
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=axisswap +order=1,-2");
+    }
+}