always_xy=True leaves a residual axisswap for Geographic3D <-> VerticalCRS pipelines (normalizeForVisualization)

[phaarnes]

Description of the issue

normalizeForVisualization() - the function behind pyproj's always_xy=True
and projinfo --normalize-axis-order is supposed to guarantee that both
ends of a coordinate operation use (lon, lat) / (E, N) axis order. For
operations between a Geographic3D CRS and a VerticalCRS (or a
CompoundCRS whose first component is a VerticalCRS), it leaves a residual
+proj=axisswap +order=2,1 in the normalized pipeline, so the output
horizontal coordinates come back lat/lon instead of lon/lat, silently,
with no error.

Root cause

A VerticalCRS has no horizontal axes, so CRS::mustAxisOrderBeSwitchedForVisualization()
(src/iso19111/crs.cpp) always returns false for it. Isolateded this is correct of cource,
since a 1D CRS has no axis order of its own to reverse.

But the actual instantiated operation (e.g. Geographic3D -> Vertical via a
geoid grid) still carries a horizontal pass-through pair internally, built
from the geographic end's own axis convention on both sides of the
pipeline (src/iso19111/operation/singleoperation.cpp, the
isGeographic3DToGravityRelatedHeight export branch). That pass-through is
symmetric and correct by itself.

CoordinateOperation::normalizeForVisualization()
(src/iso19111/operation/singleoperation.cpp) asks each end independently
whether it needs a swap: swapSource=true (geographic), swapTarget=false
(vertical, no horizontal axes to swap). It therefore prepends a reversal at
only one end, breaking the operation's internal symmetry. PROJ's
PROJ-string optimizer then cancels the leading pair (the prepended reversal
against the operation's own leading swap), but there is nothing to cancel the
operation's own trailing swap, which gets into the final pipeline.

Steps to reproduce

pyproj (forward and reverse)

from pyproj.transformer import TransformerGroup

# EPSG:4937 = ETRS89 (geographic 3D, lat/lon/h)
# EPSG:5776 = NN54 height (Norway vertical CRS)
tg = TransformerGroup("EPSG:4937", "EPSG:5776", always_xy=True, allow_ballpark=False)
t = tg.transformers[0]
print(t.to_proj4())
# +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad
#                +step +inv +proj=vgridshift +grids=no_kv_href2008a.tif +multiplier=1
#                +step +proj=unitconvert +xy_in=rad +xy_out=deg
#                +step +proj=axisswap +order=2,1         should NOT be here

print(t.transform(5.5, 58.5, 100.0))
# (58.5, 5.5, 57.087...)    lon/lat swapped on output (always_xy=True should give lon/lat)

# Reverse direction leaves the residual swap on the *input* side instead:
tg_rev = TransformerGroup("EPSG:5776", "EPSG:4937", always_xy=True, allow_ballpark=False)
print(tg_rev.transformers[0].to_proj4())
# +proj=pipeline +step +proj=axisswap +order=2,1     should NOT be here
#                +step +proj=unitconvert +xy_in=deg +xy_out=rad
#                +step +proj=vgridshift +grids=no_kv_href2008a.tif +multiplier=1
#                +step +proj=unitconvert +xy_in=rad +xy_out=deg

PROJ core (C API) - the exact function pyproj's always_xy calls

projinfo --normalize-axis-order does not reproduce this: it normalizes
the CRS objects first and searches for an operation between the already-
normalized CRSs, which for a 1D vertical target yields a symmetric pair of
swaps that cancel out. The bug only appears through the operation-level
entry point, proj_normalize_for_visualization() (PJ_OBJ.cpp), i.e. the
same call pyproj's always_xy=True makes internally:

#include <proj.h>
#include <stdio.h>

int main(void) {
    PJ_CONTEXT *ctx = proj_context_create();
    PJ *p = proj_create_crs_to_crs(ctx, "EPSG:4937", "EPSG:5776", NULL);
    PJ *norm = proj_normalize_for_visualization(ctx, p);
    printf("%s\n", proj_as_proj_string(ctx, norm, PJ_PROJ_5, NULL));
    /* +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad
       +step +inv +proj=vgridshift +grids=no_kv_href2008a.tif +multiplier=1
       +step +proj=unitconvert +xy_in=rad +xy_out=deg
       +step +proj=axisswap +order=2,1                    <- residual */
    return 0;
}

Expected behaviour

When always_xy=True / proj_normalize_for_visualization() is used, a
Geographic3D <-> VerticalCRS pipeline shouldn't have any leftover
+proj=axisswap +order=2,1 in either direction. Horizontal coordinates
should stay in (lon, lat) order all the way through, exactly as always_xy
promises. (Without it, the CRS's own declared axis order applies as usual;
this issue is only about what happens once that flag is set.)

I've got a fix for this (tested both directions, ran the full PROJ test suite
against it, no regressions). Happy to open a PR if that's useful.

Important

The same class of bug likely affects EngineeringCRS pairs with
non-trivial axis directions (e.g. EPSG:5800, axis order north/west), but
fixing that requires an axis direction flip that
Conversion::createAxisOrderReversal cannot currently express, which is a
separate, larger change. This issue and fix are scoped to the
Geographic3D <-> VerticalCRS case only.