Vertical unit sometimes not considered in transformation.

[jjimenezshaw]

Let's say we have a Geographic3D CRS, with the vertical axis in foot (in the file geog3d.txt):

$ cat geog3d.txt 
GEOGCRS["NAD83(2011) (ft)",
    DATUM["NAD83 (National Spatial Reference System 2011)",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,3],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["foot",0.3048]]

When I try to convert it to a compound system with the vertical in feet, it does no apply correctly the units conversion:

$ projinfo -s "`cat geog3d.txt`" -t EPSG:6318+8228 -o proj --spatial-test intersects --grid-check discard_missing -q
+proj=unitconvert +z_in=ft +z_out=m

If I understand this properly, it is getting the input in feet, but giving the output in meters, when I expect the output also in feet, the target CRS.
(I do not care about geoid models so far)

If source or destination is in meters and the other in feet, it works in these two scenarios:

$ projinfo -s EPSG:6319 -t EPSG:6318+8228 -o proj --spatial-test intersects --grid-check discard_missing -q
+proj=unitconvert +z_in=m +z_out=ft

$ projinfo -s "`cat geog3d.txt`" -t EPSG:6319 -o proj --spatial-test intersects --grid-check discard_missing -q
+proj=unitconvert +z_in=ft +z_out=m

However it applies the unit conversion twice in this case:

$ projinfo -s "`cat geog3d.txt`" -t EPSG:6318+5703 -o proj --spatial-test intersects --grid-check discard_missing -q
+proj=pipeline
  +step +proj=unitconvert +z_in=ft +z_out=m
  +step +proj=unitconvert +z_in=ft +z_out=m

Other wrong case I found is going to US survey foot:

$ projinfo -s "`cat geog3d.txt`" -t EPSG:6318+6360 -o proj --spatial-test intersects --grid-check discard_missing -q
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=affine +s33=1.000002000004
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=rad +z_in=ft +xy_out=deg +z_out=m

Here it's giving the output in meters, not in US survey feet. I would expect not doing +z_in=ft +z_out=m (leaving the scale for the parameter s33), or doing only +z_in=ft +z_out=us-ft without affine pipeline. More or less as in

$ projinfo -s EPSG:6319 -t EPSG:6318+6360 -o proj --spatial-test intersects --grid-check discard_missing -q
+proj=unitconvert +z_in=m +z_out=us-ft

Helper names:

• EPSG:6318 : NAD83(2011) - 2D
• EPSG:6319 : NAD83(2011) - 3D
• EPSG:8228 : NAVD88 height (ft)
• EPSG:6360 : NAVD88 height (ftUS)
• EPSG:5703 : NAVD88 height

(Note: my real use case is with a compound system with a derived vertical using the "ellipsoid" datum, but I tried to simplify the tests)

Expected behaviour

The vertical units from source and destination are consistent in the transformations.

Environment Information

• PROJ version (proj) : 9.3.0
• Operation System Information : Ubuntu 22.04

Installation method

• compiled from source