How to append z-unit conversion to Transformer pipeline?

[dmahr1]

In old versions of PyProj, transforming (x, y, z) coordinates from a 2D CRS in feet to 2D CRS in meters would automatically convert the z-axis from feet to meters as well. For example:

>>> pyproj.__version__
'1.9.6'
>>> src = pyproj.Proj(init='epsg:6578', preserve_units=True)  # state plane ft
>>> dst = pyproj.Proj(init='epsg:4326', preserve_units=True)
>>> pyproj.transform(src, dst, 3200000, 10200000, 650)
(-97.46089333100083, 30.61755306782362, 198.12039624079236)  # now in meters :)

But starting with PROJ 6.3.0 this behavior changed such that the z-axis units are not converted:

>>> pyproj.__version__
'3.0.1'
>>> transformer = pyproj.Transformer.from_crs('epsg:6578', 'epsg:4326', always_xy=True)
>>> transformer.transform(3200000, 10200000, 650)
(-97.46089333100053, 30.61755306782366, 650.0)  # still in feet :(

As Even explained in PROJ issue 2624, the z-units conversion between a pair of 2D CRSs is undefined behavior. For my purposes, I want to assume that the vertical units match the horizontal units and apply a PROJ unit conversion if they differ. I saw the discussion on GIS.SE about the to_3d() method but I don't want to do the ellipsoidal transformation, just the unit conversion.

What's the best way to insert a +proj=unitconvert to the transform pipeline? Originally I tried something like this, but it suffered from the inaccuracies associated with the PROJ4 dict:

crs = pyproj.CRS('epsg:6578')
crs_params = crs.to_dict()
if 'vunits' not in crs_params and 'units' in crs_params:
    # NOT ACCURATE, DO NOT USE
    crs = pyproj.CRS(**crs_params, vunits=crs_params['units'])

So here's an alternative approach that uses the pipelines instead of the CRS. Is this the right approach in the PyProj 2.x+ world?

transformer = pyproj.Transformer.from_crs('epsg:6578', 'epsg:4326', always_xy=True)
crs1 = pyproj.CRS('epsg:6578')
crs2 = pyproj.CRS('epsg:4326')
if (
    len(crs1.axis_info) == len(crs2.axis_info) == 2
    and crs1.axis_info[0].unit_name != crs2.axis_info[0].unit_name
):
    unit_names = ['metre' if c.is_geographic else c.axis_info[0].unit_name for c in [crs1, crs2]]
    short_names = [pyproj.database.get_units_map()[u].proj_short_name for u in unit_names]
    transformer = pyproj.Transformer.from_pipeline(
        transformer.definition + f' step proj=unitconvert z_in={short_names[0]} z_out={short_names[1]}'
    )

EDIT: I just realized this won't work when there are multiple transformations available 🤦‍♂️

[snowman2]

You could try using a compound CRS for the destination CRS if it is projected and add a VerticalCRS with units in feet: https://pyproj4.github.io/pyproj/stable/build_crs.html

For geographic CRS, you could use the to_3d() method you mentioned.

[snowman2]

Here is an example:

>>> from pyproj.crs.enums import VerticalCSAxis
>>> from pyproj.crs.coordinate_system import VerticalCS
>>> from pyproj.crs import CompoundCRS, VerticalCRS
>>> from pyproj import Transformer
>>> vcrs = VerticalCRS(
...     name="NAVD88 height",
...     datum="North American Vertical Datum 1988",
...     vertical_cs=VerticalCS(VerticalCSAxis.UP_US_FT))
>>> crs = CompoundCRS(name="epsg:6578 + navd88 height", components=["EPSG:6578", vcrs])
>>> transformer = Transformer.from_crs(crs, "EPSG:4979", always_xy=True)
>>> transformer.transform(3200000, 10200000, 650)
(-97.46089333100052, 30.61755306782416, 171.93341213002256)

[dmahr1]

Thanks for the suggestion, @snowman2! Here's what I came up with, using the epsg:{horizontal}+{vertical} syntax to avoid having to instantiate the vertical CRS manually. It assumes NAVD88 and only handles Feet and US Survey Feet. But that covers pretty much all of the cases when this system relies on the PyProj 1.x behavior that automatically converted units. Seems reasonable?

pyproj.network.set_network_enabled(True)
UNIT_TO_NAVD88_VCRS = {'foot': '+8228', 'US survey foot': '+6360'}

def get_transformer(src: str, dst: str) -> pyproj.Transformer:
    src_crs, dst_crs = pyproj.CRS(src), pyproj.CRS(dst)
    if (
        src_crs.axis_info[0].unit_name != dst_crs.axis_info[0].unit_name
        and any(crs.axis_info[0].unit_name in UNIT_TO_NAVD88_VCRS for crs in [src_crs, dst_crs])
    ):
        src += UNIT_TO_NAVD88_VCRS.get(src_crs.axis_info[0].unit_name, '+5703') if not src_crs.is_vertical else ''
        dst += UNIT_TO_NAVD88_VCRS.get(dst_crs.axis_info[0].unit_name, '+5703') if not dst_crs.is_vertical else ''
    return pyproj.Transformer.from_crs(src, dst, always_xy=True)

transformer = get_transformer('epsg:6578', 'epsg:4326')  # src will be 'epsg:6578+6360'
transformer.transform(3200000, 10200000, 650) # outputs (97.46089333100053, 30.61755306782366, 198.12039624079247)

[snowman2]

Nice 👍