Building Compound CRS with WGS84 & EGM96?

[glostis]

Disclaimer: I have read https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 and am aware of the fact that PROJ now honors the axis order of a CRS when initialized with EPSG:xxxx.

Code Sample

I think I have found an inconsistency in the way the axis order is interpreted with EPSG:4326.
Consider the following snippet:

import pyproj

ellipsoid = pyproj.CRS.from_epsg(4326)  # WGS84
geoid = pyproj.CRS.from_epsg(5773)  # EGM96
utm = pyproj.CRS.from_epsg(32610)  # UTM zone 10

lat, lon = 45, -122
for name, trf in zip(
    ["WGS84 -> EGM96", "WGS84 -> UTM"],
    [pyproj.Transformer.from_crs(ellipsoid, geoid), pyproj.Transformer.from_crs(ellipsoid, utm)],
):
    print(name, "lon, lat", trf.transform(lon, lat, 10))
    print(name, "lat, lon", trf.transform(lat, lon, 10))

It prints the following output:

WGS84 -> EGM96 lon, lat (-122.0, 45.0, 30.38683319091797)
WGS84 -> EGM96 lat, lon (45.0, -122.0, -inf)
WGS84 -> UTM lon, lat (inf, inf, inf)
WGS84 -> UTM lat, lon (578815.302916712, 4983436.768349296, 10.0)

Problem description

Since the axis order of EPSG:4326 is

Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)

I would expect to always have to pass x=lat, y=lon in a Transformer.transform call in order for everything to work correctly.

However as shown with the snippet above:

1. when transforming datums from WGS84 to EGM96, it is the call with x=lon, y=lat that gives the correct output
2. when transforming coordinates from WGS84 to UTM zone 10, as expected it is the call with x=lat, y=lon that gives the correct output.

Expected Output

I would expect to always have to pass x=lat, y=lon in Transformer.transform when my crs_from is EPSG:4326, irrespective of what my crs_to is.

Environment Information

2.6.1.post1

Installation method

pip wheel

[snowman2]

EPSG:5773 is a Vertical CRS and does not have the lat/lon component:

>>> import pyproj
>>> geoid = pyproj.CRS.from_epsg(5773)
>>> geoid
<Vertical CRS: EPSG:5773>
Name: EGM96 height
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: EGM96 geoid
- Ellipsoid: undefined
- Prime Meridian: undefined

So, you need to build the CRS you want: https://pyproj4.github.io/pyproj/stable/build_crs.html

>>> from pyproj.crs import CompoundCRS
>>> cmpd_crs = CompoundCRS(name="WGS 84 + EGM96 height", components=["EPSG:4326", "EPSG:5773"])
>>> cmpd_crs
<Compound CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: WGS GeoID
Axis Info [ellipsoidal|vertical]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- H[up]: Gravity-related height (metre)
Area of Use:
- undefined
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Sub CRS:
- WGS 84
- EGM96 height

EPSG:4326 is not 3D by nature, however EPSG:4979 is:

>>> pyproj.CRS(4979)
<Geographic 3D CRS: EPSG:4979>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: World (by country)
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

With that, your axis order should be consistent.

>>> trans = pyproj.Transformer.from_crs(cmpd_crs, "EPSG:4979")
>>> trans.transform(45, -122, 10)
(45.0, -122.0, 10.0)

[glostis]

@snowman2 thanks a lot for taking the time to write this detailed answer!

As a CRS-newbie everything PROJ-related can be a bit overwhelming...
I had not read about CompoundCRS yet, but it looks very handy! The solution you posted does exactly what I want, thanks for that 😄

Out of interest, is it expected that Transformer.transform manages to transform 3D coordinates from a 2D CRS to a Vertical CRS? Or does it work "by chance"?

---

EDIT: There seems to be a small typo in the very last line of your answer. I think it should read:

(45.0, -122.0, -10.386833190917958)

instead of

(45.0, -122.0, 10.0)

[snowman2]

Out of interest, is it expected that Transformer.transform manages to transform 3D coordinates from a 2D CRS to a Vertical CRS? Or does it work "by chance"?

It does some auto-conversions behind the scenes. But, with limited information it may or may not produce what you want. Probably where the axis-order inconsistencies come into play.

EDIT: There seems to be a small typo in the very last line of your answer. I think it should read:

Probably due to my environment missing the needed grid.

[snowman2]

Side note, here is the shorthand version for the compound CRS:

>>> cc = CRS("EPSG:4326+5773")
>>> cc
<Compound CRS: EPSG:4326+5773>
Name: WGS 84 + EGM96 height
Axis Info [ellipsoidal|vertical]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- H[up]: Gravity-related height (metre)
Area of Use:
- undefined
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Sub CRS:
- WGS 84
- EGM96 height

[glostis]

Thanks!

As far as I'm concerned the issue is resolved so feel free to close it

[elliott-imhoff]

It seems this has stopped working in the update to version 3. Running the above code prints (45.0, -122.0, 10.0) as the output (unchanged). Downgrading the pyproj version back to 2.6.1 and rerunning prints the expected (45.0, -122.0, -10.386833190917958).

[snowman2]

See: #630 (comment)

https://pyproj4.github.io/pyproj/stable/transformation_grids.html