Geoid vs ellipsoid heights for weather models

[jlmaurer]

The issue:

Weather model data is accessed through various APIs that typically provides either geopotential heights or pressure levels. These need to be converted to ellipsoidal heights rather than geoid heights.

The Good News:

After some digging, I found that RAiDER does have functionality that should convert geoid heights to ellipsoid (I think):

def _get_g_ll(lats):
    '''
    Compute the variation in gravity constant with latitude
    '''
    # TODO: verify these constants. In particular why is the reference g different from self._g0?
    return 9.80616 * (1 - 0.002637 * cosd(2 * lats) + 0.0000059 * (cosd(2 * lats))**2)


def _get_Re(lats):
    '''
    Returns: the ellipsoid as a fcn of latitude
    '''
    # TODO: verify constants, add to base class constants?
    return np.sqrt(1 / (((cosd(lats)**2) / Rmax**2) + ((sind(lats)**2) / Rmin**2)))


def _geo_to_ht(lats, hts):
    """Convert geopotential height to altitude."""
    # Convert geopotential to geometric height. This comes straight from
    # TRAIN
    # Map of g with latitude (I'm skeptical of this equation - Ray)
    g_ll = _get_g_ll(lats)
    Re = _get_Re(lats)

    # Calculate Geometric Height, h
    h = (hts * Re) / (g_ll / g0 * Re - hts)
    # from metpy
    # return (geopotential * Re) / (g0 * Re - geopotential)

    return h

The function _geo_to_ht gets called by ecmwf.py and hrrr.py to convert heights, so presumably those ones are actually being converted to ellipsoidal heights.

The Bad News

• I'm not certain these functions defined in RAiDER actually correctly convert from geoid heights to ellipsoidal heights, or which ellipsoid they convert to, or why metpy has a different function than TRAIN
• Not all the weather models use these functions. ecmwf.py and hrrr.py use it, but GMAO, NCMR, MERRA-2 do not. It does look like all the ECMWF models call the function
• It's not 100% clear to me which geoid the models use. In this issue I reference the ECMWF page that states geoidal heights are in EGM96, but not sure if this is true for other models, does it matter, etc.

[jlmaurer]

Pinging @dbekaert @piyushrpt @bbuzz31 .
@dbekaert I was wrong...at least ECMWF models and HRRR seem to be getting converted to ellipsoidal heights. However, the other models are not, and I'm not sure which ellipsoid this one is.

[piyushrpt]

1. Weather models typically do not use an ellipsoid model for the earth. They usually use a spherical model - you can look at the proj strings and see that a radius is usually provided and no eccentricity. Geopotential heights are defined on top of this sphere - the equivalent would be ellipsoidal heights.
2. When ppl try to use weather models with other data, they tend to convert cartesian coordinates to spherical lat/lon and then assume that lat/lon represents points on the ellipsoid. So, the whole handling of map projections is a bit of a dance. Assumption is that these usually make very little difference due to low resolution of weather models. ECMWF / HRRR - we assign EPSG:4326 I guess but if you were to look up ECMWF documentation that might tell you what spherical radius was used - this is embedded in some code in the grib file.
3. All geometric computations that have closed form expressions use ellipsoidal heights.

[piyushrpt]

https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation#ERA5:datadocumentation-Earthmodel

[jlmaurer]

Thanks @piyushrpt ! Just to document specifically the conclusions for our code:

1. Global weather models typically use a sphere for their calculations, not an ellipsoid, but we can assume that spherical lat/lons are equal to geodetic lat/lons and not be too far off (longitude should in general be fine; latitudes will have a latitude-dependent error that is largest +/- 45 degrees from the equator).
2. In RAiDER, geopotential heights get converted to heights above the WGS84 ellipsoid for both ECMWF models and HRRR. This done by using a latitude-dependent g and Earth radius. So, weather model heights in the processed weather model are in WGS84 heights, or at least as close as you can get with the approximations made.
3. We need to double-check GMAO and MERRA-2 to see what they use.
4. I've updated the docstrings for the relevant functions to (hopefully) clarify what is what and also link to external sites for reference.

[piyushrpt]

Essentially (lat,lon,hgt) on sphere is assumed to be the same as (lat, lon, hgt) on the ellipsoid. With these other models, try tranforming (Xmap, Ymap, ht) to (Lon, Lat, ht) and if the ht remains the same - the path being used for coordinate transform does not worry about difference between sphere and ellipsoid.

To be super precise you may want:

ECEF (4978) -> LatLon (4979) - take the transformed array as is and feed it to LatLon (sphere) -> ModelCRS .

We are currently using: ECEF (4978) -> ModelCRS directly.

This is all from an old blog post about WRF: https://fabienmaussion.info/2018/01/06/wrf-projection/

[piyushrpt]

It is this specific line that needs to be replaced with 2 calls:

RAiDER/tools/RAiDER/delay.py

Line 389 in 2f28cb5

pts = ecef_to_model.transform(

ecef2lla -> pass lla array to -> llasphere_to_model

you will need to build llasphere using the radius value of the model and proj strings like +proj=latlon +R=value_from_model +no_defs