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

New functions for calculation of geopotential from height and vice versa #678

Merged
merged 9 commits into from Dec 18, 2017
Merged

New functions for calculation of geopotential from height and vice versa #678

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
ad7c444
Calculation of geopotential and geopot. height
Dec 8, 2017
696bbcb
New function for conversion from geopot to height
Dec 11, 2017
145ed1e
Update and style changes to geopotential functions
Dec 12, 2017
6e81070
Add tests for geopotential functions and fix typo
Dec 12, 2017
c9b3df9
Replace single with double quotes in docstrings
Dec 12, 2017
1d5ae99
Fix alphabetical import order
Dec 13, 2017
237e8e6
Fix lint error
Dec 13, 2017
6f90805
Fix the docstring test failure
Dec 13, 2017
51a9a6d
Parenthesis on line continuation in "constants"
Dec 18, 2017
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
83 changes: 82 additions & 1 deletion metpy/calc/basic.py
View file
Open in desktop
Expand Up @@ -16,7 +16,7 @@

import numpy as np

from ..constants import g, omega, Rd
from ..constants import G, g, me, omega, Rd, Re
from ..package_tools import Exporter
from ..units import atleast_1d, check_units, masked_array, units

Expand Down Expand Up @@ -257,6 +257,87 @@ def pressure_to_height_std(pressure):
return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(Rd * gamma / g))


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

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

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

Examples
--------
>>> from metpy.constants import g, G, me, Re
>>> import metpy.calc
>>> from metpy.units import units
>>> height = np.linspace(0,10000, num = 11) * units.m
>>> geopot = 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')>

Notes
-----
Derived from definition of geopotential in [Hobbs2006]_ pg.14 Eq.1.8.

"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should be a blank line at the end of the docstring (inside the quotes).

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added the blank lines in the docstring(s).

# Calculate geopotential
geopot = G * me * ((1 / Re) - (1 / (Re + height)))

return geopot


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

Parameters
----------
geopotential : `pint.Quantity`
Geopotential (array_like)

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

Examples
--------
>>> from metpy.constants import g, G, me, Re
>>> import metpy.calc
>>> from metpy.units import units
>>> height = np.linspace(0,10000, num = 11) * units.m
>>> geopot = 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')>
>>> height = metpy.calc.geopotential_to_height(geopot)
>>> height
<Quantity([ 0. 1000. 2000. 3000. 4000. 5000. 6000. 7000.
8000. 9000. 10000.], 'meter')>

Notes
-----
Derived from definition of geopotential in [Hobbs2006]_ pg.14 Eq.1.8.

"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Blank line.

# Calculate geopotential
height = (((1 / Re) - (geopot / (G * me))) ** -1) - Re

return height


@exporter.export
@check_units('[length]')
def height_to_pressure_std(height):
Expand Down
22 changes: 19 additions & 3 deletions metpy/calc/tests/test_basic.py
View file
Open in desktop
Expand Up @@ -7,9 +7,9 @@
import pytest

from metpy.calc import (add_height_to_pressure, add_pressure_to_height, coriolis_parameter,
get_wind_components, get_wind_dir, get_wind_speed, heat_index,
height_to_pressure_std, pressure_to_height_std, sigma_to_pressure,
windchill)
geopotential_to_height, get_wind_components, get_wind_dir, get_wind_speed,
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

E501 line too long (98 > 95 characters)

heat_index, height_to_geopotential, height_to_pressure_std,
pressure_to_height_std, sigma_to_pressure, windchill)
from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_array_equal
from metpy.units import units

Expand Down Expand Up @@ -194,6 +194,22 @@ def test_heat_index_ratio():
hi = heat_index(temp, rh)
assert_almost_equal(hi.to('degC'), units.Quantity([50.3405, np.nan], units.degC), 4)


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)


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'))
height = geopotential_to_height(geopotential)
assert_array_almost_equal(height, units.Quantity([0, 1000, 2000, 3000], units.m), 0)

# class TestIrrad(object):
# def test_basic(self):
# 'Test the basic solar irradiance calculation.'
Expand Down
28 changes: 17 additions & 11 deletions metpy/constants.py
View file
Open in desktop
Expand Up @@ -6,17 +6,19 @@

Earth
-----
======================== =============== =========== ========================== =======================================
Name Symbol Short Name Units Description
------------------------ --------------- ----------- -------------------------- ---------------------------------------
earth_avg_radius :math:`R_e` Re :math:`\text{m}` Avg. radius of the Earth
earth_gravity :math:`g` g :math:`\text{m s}^{-2}` Avg. gravity acceleration on Earth
earth_avg_angular_vel :math:`\Omega` omega :math:`\text{rad s}^{-1}` Avg. angular velocity of Earth
earth_sfc_avg_dist_sun :math:`d` d :math:`\text{m}` Avg. distance of the Earth from the Sun
earth_solar_irradiance :math:`S` S :math:`\text{W m}^{-2}` Avg. solar irradiance of Earth
earth_max_declination :math:`\delta` delta :math:`\text{degrees}` Max. solar declination angle of Earth
earth_orbit_eccentricity :math:`e` :math:`\text{None}` Avg. eccentricity of Earth's orbit
======================== =============== =========== ========================== =======================================
======================== =============== =========== ======================================= ========================================
Name Symbol Short Name Units Description
------------------------ --------------- ----------- --------------------------------------- ----------------------------------------
earth_avg_radius :math:`R_e` Re :math:`\text{m}` Avg. radius of the Earth
earth_gravity :math:`g` g :math:`\text{m s}^{-2}` Avg. gravity acceleration on Earth
gravitational_constant :math:`G` G :math:`\text{m}^{3} {kg}^{-1} {s}^{-2}` Gravitational constant
earth_avg_angular_vel :math:`\Omega` omega :math:`\text{rad s}^{-1}` Avg. angular velocity of Earth
earth_sfc_avg_dist_sun :math:`d` d :math:`\text{m}` Avg. distance of the Earth from the Sun
earth_solar_irradiance :math:`S` S :math:`\text{W m}^{-2}` Avg. solar irradiance of Earth
earth_max_declination :math:`\delta` delta :math:`\text{degrees}` Max. solar declination angle of Earth
earth_orbit_eccentricity :math:`e` :math:`\text{None}` Avg. eccentricity of Earth's orbit
earth_mass :math:`m_e` me :math:`\text{kg}` Total mass of the Earth (approx)
======================== =============== =========== ======================================= ========================================

Water
-----
Expand Down Expand Up @@ -70,11 +72,15 @@
earth_gravity = g = units.Quantity(1.0, units.gravity).to('m / s^2')
# Taken from GEMPAK constants
Re = earth_avg_radius = 6.3712e6 * units.m
G = gravitational_constant = units.Quantity(1, units.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you remove the \ and just wrap the whole thing in parenthesis?

newtonian_constant_of_gravitation) \
.to('m^3 / kg / s^2')
omega = earth_avg_angular_vel = 2 * units.pi / units.sidereal_day
d = earth_sfc_avg_dist_sun = 1.496e11 * units.m
S = earth_solar_irradiance = units.Quantity(1.368e3, 'W / m^2')
delta = earth_max_declination = 23.45 * units.deg
earth_orbit_eccentricity = 0.0167
earth_mass = me = 5.9722e24 * units.kg

# molar gas constant
R = units.Quantity(1.0, units.R).to('J / K / mol')
Expand Down