createOperations(): fix GeogCRS to GeogCRS when there is a axis reversal and vertical unit change (relates to https://github.com/OSGeo/gdal/issues/1561)

2 files changed, 64 insertions, 28 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 8362b3a2..3b480a92 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11302,12 +11302,38 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
         geogSrc->datum()->_isEquivalentTo(
             geogDst->datum().get(), util::IComparable::Criterion::EQUIVALENT);
 
+    // Do the CRS differ by their axis order ?
+    bool axisReversal2D = false;
+    bool axisReversal3D = false;
+    if (!srcCS->_isEquivalentTo(dstCS.get(),
+                                util::IComparable::Criterion::EQUIVALENT)) {
+        auto srcOrder = srcCS->axisOrder();
+        auto dstOrder = dstCS->axisOrder();
+        if (((srcOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST ||
+              srcOrder == cs::EllipsoidalCS::AxisOrder::
+                              LAT_NORTH_LONG_EAST_HEIGHT_UP) &&
+             (dstOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH ||
+              dstOrder == cs::EllipsoidalCS::AxisOrder::
+                              LONG_EAST_LAT_NORTH_HEIGHT_UP)) ||
+            ((srcOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH ||
+              srcOrder == cs::EllipsoidalCS::AxisOrder::
+                              LONG_EAST_LAT_NORTH_HEIGHT_UP) &&
+             (dstOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST ||
+              dstOrder == cs::EllipsoidalCS::AxisOrder::
+                              LAT_NORTH_LONG_EAST_HEIGHT_UP))) {
+            if (srcAxisList.size() == 3 || dstAxisList.size() == 3)
+                axisReversal3D = true;
+            else
+                axisReversal2D = true;
+        }
+    }
+
     // Do they differ by vertical units ?
     if (vconvSrc != vconvDst &&
         geogSrc->ellipsoid()->_isEquivalentTo(
             geogDst->ellipsoid().get(),
             util::IComparable::Criterion::EQUIVALENT)) {
-        if (offset_pm.value() == 0) {
+        if (offset_pm.value() == 0 && !axisReversal2D && !axisReversal3D) {
             // If only by vertical units, use a Change of Vertical
             // Unit
             // transformation
@@ -11329,33 +11355,11 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog(
     }
 
     // Do the CRS differ only by their axis order ?
-    if (sameDatum &&
-        !srcCS->_isEquivalentTo(dstCS.get(),
-                                util::IComparable::Criterion::EQUIVALENT)) {
-        auto srcOrder = srcCS->axisOrder();
-        auto dstOrder = dstCS->axisOrder();
-        if ((srcOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST &&
-             dstOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH) ||
-            (srcOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH &&
-             dstOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST)) {
-            auto conv = Conversion::createAxisOrderReversal(false);
-            conv->setCRSs(sourceCRS, targetCRS, nullptr);
-            res.emplace_back(conv);
-            return res;
-        }
-        if ((srcOrder ==
-                 cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST_HEIGHT_UP &&
-             dstOrder ==
-                 cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH_HEIGHT_UP) ||
-            (srcOrder ==
-                 cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH_HEIGHT_UP &&
-             dstOrder ==
-                 cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST_HEIGHT_UP)) {
-            auto conv = Conversion::createAxisOrderReversal(true);
-            conv->setCRSs(sourceCRS, targetCRS, nullptr);
-            res.emplace_back(conv);
-            return res;
-        }
+    if (sameDatum && (axisReversal2D || axisReversal3D)) {
+        auto conv = Conversion::createAxisOrderReversal(axisReversal3D);
+        conv->setCRSs(sourceCRS, targetCRS, nullptr);
+        res.emplace_back(conv);
+        return res;
     }
 
     std::vector<CoordinateOperationNNPtr> steps;
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 2ba76172..64112a82 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4759,6 +4759,38 @@ TEST(operation, geogCRS_to_geogCRS_3D) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_3D_lat_long_non_metre_to_geogCRS_longlat) {
+
+    auto wkt = "GEOGCRS[\"my CRS\",\n"
+               "    DATUM[\"World Geodetic System 1984\",\n"
+               "        ELLIPSOID[\"WGS 84\",6378137,298.257223563],\n"
+               "        ID[\"EPSG\",6326]],\n"
+               "    CS[ellipsoidal,3],\n"
+               "        AXIS[\"latitude\",north,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+               "        AXIS[\"longitude\",east,\n"
+               "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+               "        AXIS[\"ellipsoidal height\",up,\n"
+               "            LENGTHUNIT[\"my_vunit\",0.3]]]";
+    auto srcCRS_obj = WKTParser().createFromWKT(wkt);
+    auto srcCRS = nn_dynamic_pointer_cast<CRS>(srcCRS_obj);
+    ASSERT_TRUE(srcCRS != nullptr);
+
+    auto dstCRS_obj = PROJStringParser().createFromPROJString(
+        "+proj=longlat +datum=WGS84 +type=crs");
+    auto dstCRS = nn_dynamic_pointer_cast<CRS>(dstCRS_obj);
+    ASSERT_TRUE(dstCRS != nullptr);
+
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(srcCRS), NN_CHECK_ASSERT(dstCRS));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+              "+proj=unitconvert +z_in=0.3 +z_out=m");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, geogCRS_without_id_to_geogCRS_3D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");