Getting geoid height

[NilsKeunecke]

I want to extract geoid height at a location (lat, lon) to get from EGM96 elevation to WGS84 ellipsoid elevation.

Essentially solve for N in the equation here.

First I define by standard WGS84 ellipsoid CRS:

from pyproj import Transformer, CRS
from pyproj.crs import CompoundCRS

wgs84_ellipsoid_crs = CRS("EPSG: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.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

And secondly, I define a combination of a 2D WGS84 with EGM96 elevation:

wgs84_no_altitude_crs = CRS("EPSG:4326")# 2D WGS84

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich


egm96 = CRS("EPSG:5773") # Vertical component geoid altitude

 <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

And then I finally combine by 2D CRS with my Vertical CRS:

wgs84_egm96_crs = CompoundCRS(name="wgs84_egm96", components=["EPSG:4326", "EPSG:5733"])

Even though the egm96 and wgs84_no_altitude share the same area of use, the CompoundCRS seems to have no area of use:

<Compound CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: wgs84_egm96_crs
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 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Sub CRS:
- WGS 84
- DNN height

I saw that one can inspect if there is missing information for the transformation and download it.:

from pyproj.transformer import TransformerGroup
tg = TransformerGroup(wgs84_crs, wgs84_egm96_crs_crs)
tg.download_grids(verbose=True)
<TransformerGroup: best_available=True>
- transformers: 2
- unavailable_operations: 0

This indicates that all should be good to go. But when I try to use the transformation, it does nothing:

lat, lon, alt = 48.85813557146666, 2.338315412849036, 33.8
trans = Transformer.from_crs(wgs84_crs,wgs84_egm96_crs_crs)
trans.transform(lat, lon, alt, inplace=True) # Expecting an altitude of 44.53

Questions:
Why is the area of use undefined for wgs84_egm96_crs?
How do I combine the two CRS properly? (How do I combine a Geographic 2D CRS with a Vertical CRS?)

[j-carson]

EPSG 9707 has been working for me: https://epsg.io/9707
Note that when using download_grids, it may not work the first time #1474