fix residual axisswap for Geographic3D <-> VerticalCRS pipelines

include/proj/crs.hpp

@@ -143,6 +143,8 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage,
 
     PROJ_INTERNAL bool mustAxisOrderBeSwitchedForVisualization() const;
 
+    PROJ_INTERNAL bool hasNoHorizontalAxes() const;
+
     PROJ_INTERNAL CRSNNPtr applyAxisOrderReversal(const char *nameSuffix) const;
 
     PROJ_FOR_TEST CRSNNPtr normalizeForVisualization() const;

src/iso19111/crs.cpp

@@ -946,6 +946,31 @@ bool CRS::mustAxisOrderBeSwitchedForVisualization() const {
     return false;
 }
 
+// ---------------------------------------------------------------------------
+
+/** Return whether this CRS declares no horizontal axes.
+ *
+ * Coordinate operations involving such a CRS just pass horizontal coordinates
+ * through, so their axis order is dictated by the other end of the operation
+ * rather than by this CRS.
+ */
+bool CRS::hasNoHorizontalAxes() const {
+
+    if (dynamic_cast<const VerticalCRS *>(this)) {
+        return true;
+    }
+
+    if (const CompoundCRS *compoundCRS =
+            dynamic_cast<const CompoundCRS *>(this)) {
+        const auto &comps = compoundCRS->componentReferenceSystems();
+        if (!comps.empty()) {
+            return comps[0]->hasNoHorizontalAxes();
+        }
+    }
+
+    return false;
+}
+
 //! @endcond
 
 // ---------------------------------------------------------------------------

src/iso19111/operation/singleoperation.cpp

@@ -422,10 +422,20 @@ CoordinateOperation::normalizeForVisualization() const {
         throw util::UnsupportedOperationException(
             "Cannot retrieve source or target CRS");
     }
-    const bool swapSource =
-        l_sourceCRS->mustAxisOrderBeSwitchedForVisualization();
-    const bool swapTarget =
-        l_targetCRS->mustAxisOrderBeSwitchedForVisualization();
+    bool swapSource = l_sourceCRS->mustAxisOrderBeSwitchedForVisualization();
+    bool swapTarget = l_targetCRS->mustAxisOrderBeSwitchedForVisualization();
+
+    // A CRS without horizontal axes lets horizontal coordinates pass through
+    // in the axis order of the other end, so both ends must be swapped
+    // together to avoid leaving a residual axis swap in the pipeline.
+    const bool sourceHasNoHorizontalAxes = l_sourceCRS->hasNoHorizontalAxes();
+    const bool targetHasNoHorizontalAxes = l_targetCRS->hasNoHorizontalAxes();
+    if (targetHasNoHorizontalAxes && !sourceHasNoHorizontalAxes) {
+        swapTarget = swapSource;
+    } else if (sourceHasNoHorizontalAxes && !targetHasNoHorizontalAxes) {
+        swapSource = swapTarget;
+    }
+
     auto l_this = NN_NO_CHECK(std::dynamic_pointer_cast<CoordinateOperation>(
         shared_from_this().as_nullable()));
     if (!swapSource && !swapTarget) {

test/unit/test_operation.cpp

@@ -5716,6 +5716,65 @@ TEST(operation, normalizeForVisualization) {
             "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
     }
 
+    // Source(geographic 3D) and target(vertical): the vertical CRS has no
+    // horizontal axes, so the horizontal coordinates it passes through must
+    // follow the source axis order. No residual axisswap must remain.
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        // ETRS89 (geographic 3D)
+        auto src = authFactory->createCoordinateReferenceSystem("4937");
+        // NN54 height
+        auto dst = authFactory->createCoordinateReferenceSystem("5776");
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            src, dst, ctxt);
+        ASSERT_GE(list.size(), 1U);
+        auto op = list[0];
+        // Without normalization the operation carries its own axis swaps.
+        EXPECT_TRUE(
+            op->exportToPROJString(PROJStringFormatter::create(
+                                       PROJStringFormatter::Convention::PROJ_5,
+                                       authFactory->databaseContext())
+                                       .get())
+                .find("+proj=axisswap +order=2,1") != std::string::npos);
+        auto opNormalized = op->normalizeForVisualization();
+        EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get()));
+        const auto projString = opNormalized->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
+                                        authFactory->databaseContext())
+                .get());
+        EXPECT_TRUE(projString.find("+proj=axisswap +order=2,1") ==
+                    std::string::npos)
+            << projString;
+        EXPECT_TRUE(projString.find("+proj=vgridshift") != std::string::npos)
+            << projString;
+    }
+
+    // Reverse of above: source(vertical) and target(geographic 3D)
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        // NN54 height
+        auto src = authFactory->createCoordinateReferenceSystem("5776");
+        // ETRS89 (geographic 3D)
+        auto dst = authFactory->createCoordinateReferenceSystem("4937");
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            src, dst, ctxt);
+        ASSERT_GE(list.size(), 1U);
+        auto op = list[0];
+        auto opNormalized = op->normalizeForVisualization();
+        EXPECT_FALSE(opNormalized->_isEquivalentTo(op.get()));
+        const auto projString = opNormalized->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5,
+                                        authFactory->databaseContext())
+                .get());
+        EXPECT_TRUE(projString.find("+proj=axisswap +order=2,1") ==
+                    std::string::npos)
+            << projString;
+        EXPECT_TRUE(projString.find("+proj=vgridshift") != std::string::npos)
+            << projString;
+    }
+
     // Source(boundCRS) and target(geographic) must be inverted
     {
         auto src = BoundCRS::createFromTOWGS84(