Add geopotential to/from height calculation

[jrleeman]

Add a simple conversion for geopotential height.

[stefan-hofer]

Combining the altitude dependency of gravity (https://en.wikipedia.org/wiki/Gravity_of_Earth#Altitude) and the equation for geopotential height (http://www.atms.unca.edu/chennon/notes/atms310/geopotential.pdf)? This seems to me a fairly straightforward problem and I'd like to help (my first pr/Issue). Would this function be added to https://github.com/Unidata/MetPy/blob/master/metpy/calc/basic.py or somewhere completely different?

[jrleeman]

Hi @shofer16450 - Great! We're always happy to welcome new contributions! Yep, that's the idea. Given that we'll want to go both to/from Geopotential height, I think I'd break it down as follows:

• Function for gravity at an altitude. I think I'd add the mass of the earth to constants and solve the gravitational equation instead of using the approximate form. Though we could also just calculate the free air anomaly knowing .3086 mgal/m decrease in gravity. That's faster, but doesn't give us the extra constants, etc. Thoughts @dopplershift ?
• Function to go to geopotential
• Function to go from geopotential

I think basic is a perfectly fine place for that.

[stefan-hofer]

I'm trying to implement the geopotential functionality now. Having done the calculation with the functions for the height dependency of gravity I think for meteorological purposes we could even assume g to be constant throughout the atmosphere, given that we're almost exclusively dealing with tropospheric (<= 15 km) problems. The error we would make with this assumption is minimal.

Given all the constants and the function for gravitational acceleration, the difference between the gravity at 10 km and at sea level for example is ~0.3%.

g_top = 9.78542882  # This is g at 10 km
g_surf = 9.81924435 # This is g at the surface (not defined 'standard gravity')

(1 - g_top/g_bottom) * 100
0.3443801660766299 # Difference in percent

It would definitely be quicker.

[stefan-hofer]

On the other hand, the code runs quite quick and the integral for geopotential is quite easy to solve. This would be my conversion from height to geopotential and geopotential height conversion. What do you think? For this solution I have added G (Gravitational constant) and the mass of the earth "me" to the constants and used the solution for the integral which is for example given here: https://en.wikipedia.org/wiki/Geopotential

@exporter.export
@check_units('[length]')
def height_to_geopotential(height):
    r"""Computes geopotential and geopotential height for a given height in meters

    Parameters
    ----------
    height : `pint.Quantity`
        Height above sea level (array_like)

    Returns
    -------
    `pint.Quantity`
        The corresponding geopotential value(s)
    `pint.Qantity`
        The corresponding geopotential height value(s)

    Examples
    --------
    >>> from metpy.constants import g,Re,G,me
    >>> import metpy.calc
    >>> from metpy.units import units
    >>> height = np.linspace(0,10000, num = 11) * units.m
    >>> geopot,geopot_height = metpy.calc.height_to_geopotential(height)
    >>> geopot
    <Quantity([     0.           9817.70342881  19632.32592389  29443.86893527
    39252.33391207  49057.72230251  58860.0355539   68659.27511263
    78455.4424242   88248.53893319  98038.56608327],
    'meter ** 2 / second ** 2')>

    >>> geopot_height
    <Quantity([    0.          1001.12713606  2001.94010431  3002.43905261
    4002.62412874 5002.49548036  6002.05325508  7001.29760037  8000.22866363
    8998.84659218  9997.15153322], 'meter')>

    """
    # Calculate geopotential
    geopot = G * me * ((1 / Re) - (1 / (Re + height)))
    # Normalize geopotential by standard gravitational acceleration
    geopot_height = geopot / g

    return geopot, geopot_height

[dopplershift]

That looks great. The only thing that gives me pause is that I lean towards only returning geopotential from the function. If we need a function for geopotential -> geopotential height (divide by g), I think that should be its own function, since that's a completely separate calculation.

[stefan-hofer]

That makes sense to me. I'll change the function so it will only return geopotential and will also write the function to get from geopotential to height.

[jrleeman]

This was closed #678