Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Change geopotential <-> height formulas to use standard gravity #1174

Merged
merged 1 commit into from Dec 24, 2019
Merged

Change geopotential <-> height formulas to use standard gravity #1174

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
81 changes: 55 additions & 26 deletions src/metpy/calc/basic.py
View file
Open in desktop
Expand Up @@ -420,16 +420,16 @@ def pressure_to_height_std(pressure):
@preprocess_xarray
@check_units('[length]')
def height_to_geopotential(height):
r"""Compute geopotential for a given height.
r"""Compute geopotential for a given height above sea level.

Calculates the geopotential from height using the following formula, which is derived from
the definition of geopotential as given in [Hobbs2006]_ Pg. 69 Eq 3.21:
Calculates the geopotential from height above mean sea level using the following formula,
which is derived from the definition of geopotential as given in [Hobbs2006]_ Pg. 69 Eq
3.21, along with an approximation for variation of gravity with altitude:

.. math:: \Phi = G m_e \left( \frac{1}{R_e} - \frac{1}{R_e + z}\right)
.. math:: \Phi = \frac{g R_e z}{R_e + z}

(where :math:`\Phi` is geopotential, :math:`z` is height, :math:`R_e` is average Earth
radius, :math:`G` is the (universal) gravitational constant, and :math:`m_e` is the
approximate mass of Earth.)
radius, and :math:`g` is standard gravity.)

Parameters
----------
Expand All @@ -448,29 +448,44 @@ def height_to_geopotential(height):
>>> height = np.linspace(0, 10000, num=11) * units.m
>>> geopot = metpy.calc.height_to_geopotential(height)
>>> geopot
<Quantity([ 0. 9817.46806283 19631.85526579 29443.16305887
39251.39289118 49056.54621087 58858.62446524 68657.62910064
78453.56156252 88246.42329544 98036.21574305], 'meter ** 2 / second ** 2')>
<Quantity([ 0. 9805.11102602 19607.14506998 29406.10358006
39201.98800351 48994.79978671 58784.54037509 68571.21121319
78354.81374467 88135.34941224 97912.81965774], 'meter ** 2 / second ** 2')>

Notes
-----
This calculation approximates :math:`g(z)` as

.. math:: g(z) = g_0 \left( \frac{R_e}{R_e + z} \right)^2

where :math:`g_0` is standard gravity. It thereby accounts for the average effects of
centrifugal force on apparent gravity, but neglects latitudinal variations due to
centrifugal force and Earth's eccentricity.

(Prior to MetPy v0.11, this formula instead calculated :math:`g(z)` from Newton's Law of
Gravitation assuming a spherical Earth and no centrifugal force effects.)

See Also
--------
geopotential_to_height

"""
# Direct implementation of formula from Hobbs yields poor numerical results (see
# gh-1075), so was replaced with algebraic equivalent.
return (mpconsts.G * mpconsts.me / mpconsts.Re) * (height / (mpconsts.Re + height))
return (mpconsts.g * mpconsts.Re * height) / (mpconsts.Re + height)


@exporter.export
@preprocess_xarray
def geopotential_to_height(geopot):
r"""Compute height from a given geopotential.
r"""Compute height above sea level from a given geopotential.

Calculates the height from geopotential using the following formula, which is derived from
the definition of geopotential as given in [Hobbs2006]_ Pg. 69 Eq 3.21:
Calculates the height above mean sea level from geopotential using the following formula,
which is derived from the definition of geopotential as given in [Hobbs2006]_ Pg. 69 Eq
3.21, along with an approximation for variation of gravity with altitude:

.. math:: z = \frac{1}{\frac{1}{R_e} - \frac{\Phi}{G m_e}} - R_e
.. math:: z = \frac{\Phi R_e}{gR_e - \Phi}

(where :math:`\Phi` is geopotential, :math:`z` is height, :math:`R_e` is average Earth
radius, :math:`G` is the (universal) gravitational constant, and :math:`m_e` is the
approximate mass of Earth.)
radius, and :math:`g` is standard gravity.)

Parameters
----------
Expand All @@ -480,7 +495,7 @@ def geopotential_to_height(geopot):
Returns
-------
`pint.Quantity`
The corresponding height value(s)
The corresponding value(s) of height above sea level

Examples
--------
Expand All @@ -489,19 +504,33 @@ def geopotential_to_height(geopot):
>>> height = np.linspace(0, 10000, num=11) * units.m
>>> geopot = metpy.calc.height_to_geopotential(height)
>>> geopot
<Quantity([ 0. 9817.46806283 19631.85526579 29443.16305887
39251.39289118 49056.54621087 58858.62446524 68657.62910064
78453.56156252 88246.42329544 98036.21574305], 'meter ** 2 / second ** 2')>
<Quantity([ 0. 9805.11102602 19607.14506998 29406.10358006
39201.98800351 48994.79978671 58784.54037509 68571.21121319
78354.81374467 88135.34941224 97912.81965774], 'meter ** 2 / second ** 2')>
>>> height = metpy.calc.geopotential_to_height(geopot)
>>> height
<Quantity([ 0. 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000.
9000. 10000.], 'meter')>

Notes
-----
This calculation approximates :math:`g(z)` as

.. math:: g(z) = g_0 \left( \frac{R_e}{R_e + z} \right)^2

where :math:`g_0` is standard gravity. It thereby accounts for the average effects of
centrifugal force on apparent gravity, but neglects latitudinal variations due to
centrifugal force and Earth's eccentricity.

(Prior to MetPy v0.11, this formula instead calculated :math:`g(z)` from Newton's Law of
Gravitation assuming a spherical Earth and no centrifugal force effects.)

See Also
--------
height_to_geopotential

"""
# Direct implementation of formula from Hobbs yields poor numerical results (see
# gh-1075), so was replaced with algebraic equivalent.
scaled = geopot * mpconsts.Re
return scaled * mpconsts.Re / (mpconsts.G * mpconsts.me - scaled)
return (geopot * mpconsts.Re) / (mpconsts.g * mpconsts.Re - geopot)


@exporter.export
Expand Down
20 changes: 10 additions & 10 deletions tests/calc/test_basic.py
View file
Open in desktop
Expand Up @@ -265,24 +265,24 @@ def test_height_to_geopotential():
"""Test conversion from height to geopotential."""
height = units.Quantity([0, 1000, 2000, 3000], units.m)
geopot = height_to_geopotential(height)
assert_array_almost_equal(geopot, units.Quantity([0., 9817, 19632,
29443], units('m**2 / second**2')), 0)
assert_array_almost_equal(geopot, units.Quantity([0., 9805, 19607,
29406], units('m**2 / second**2')), 0)


# See #1075 regarding previous destructive cancellation in floating point
def test_height_to_geopotential_32bit():
"""Test conversion to geopotential with 32-bit values."""
heights = np.linspace(20597, 20598, 11, dtype=np.float32) * units.m
truth = np.array([201590.422, 201591.391, 201592.375, 201593.344,
201594.312, 201595.297, 201596.266, 201597.250,
201598.219, 201599.203, 201600.172], dtype=np.float32) * units('J/kg')
truth = np.array([201336.67, 201337.66, 201338.62, 201339.61, 201340.58, 201341.56,
201342.53, 201343.52, 201344.48, 201345.44, 201346.42],
dtype=np.float32) * units('J/kg')
assert_almost_equal(height_to_geopotential(heights), truth, 2)


def test_geopotential_to_height():
"""Test conversion from geopotential to height."""
geopotential = units.Quantity([0, 9817.70342881, 19632.32592389,
29443.86893527], units('m**2 / second**2'))
geopotential = units.Quantity([0., 9805.11102602, 19607.14506998, 29406.10358006],
units('m**2 / second**2'))
height = geopotential_to_height(geopotential)
assert_array_almost_equal(height, units.Quantity([0, 1000, 2000, 3000], units.m), 0)

Expand All @@ -291,9 +291,9 @@ def test_geopotential_to_height():
def test_geopotential_to_height_32bit():
"""Test conversion from geopotential to height with 32-bit values."""
geopot = np.arange(201590, 201600, dtype=np.float32) * units('J/kg')
truth = np.array([20596.957, 20597.059, 20597.162, 20597.266,
20597.365, 20597.469, 20597.570, 20597.674,
20597.777, 20597.881], dtype=np.float32) * units.m
truth = np.array([20623.000, 20623.102, 20623.203, 20623.307, 20623.408,
20623.512, 20623.615, 20623.717, 20623.820, 20623.924],
dtype=np.float32) * units.m
assert_almost_equal(geopotential_to_height(geopot), truth, 2)


Expand Down