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 3 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
135 changes: 108 additions & 27 deletions metpy/calc/basic.py
View file
Open in desktop
@@ -1,22 +1,22 @@
# Copyright (c) 2008,2015,2016 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains a collection of basic calculations.
'''Contains a collection of basic calculations.

These include:

* wind components
* heat index
* windchill
"""
'''

from __future__ import division

import warnings

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 All @@ -25,7 +25,7 @@

@exporter.export
def get_wind_speed(u, v):
r"""Compute the wind speed from u and v-components.
r'''Compute the wind speed from u and v-components.

Parameters
----------
Expand All @@ -43,14 +43,14 @@ def get_wind_speed(u, v):
--------
get_wind_components

"""
'''
speed = np.sqrt(u * u + v * v)
return speed


@exporter.export
def get_wind_dir(u, v):
r"""Compute the wind direction from u and v-components.
r'''Compute the wind direction from u and v-components.

Parameters
----------
Expand All @@ -69,7 +69,7 @@ def get_wind_dir(u, v):
--------
get_wind_components

"""
'''
wdir = 90. * units.deg - np.arctan2(-v, -u)
origshape = wdir.shape
wdir = atleast_1d(wdir)
Expand All @@ -79,7 +79,7 @@ def get_wind_dir(u, v):

@exporter.export
def get_wind_components(speed, wdir):
r"""Calculate the U, V wind vector components from the speed and direction.
r'''Calculate the U, V wind vector components from the speed and direction.

Parameters
----------
Expand Down Expand Up @@ -107,7 +107,7 @@ def get_wind_components(speed, wdir):
(<Quantity(7.071067811865475, 'meter / second')>,
<Quantity(7.071067811865477, 'meter / second')>)

"""
'''
wdir = _check_radians(wdir, max_radians=4 * np.pi)
u = -speed * np.sin(wdir)
v = -speed * np.cos(wdir)
Expand All @@ -117,7 +117,7 @@ def get_wind_components(speed, wdir):
@exporter.export
@check_units(temperature='[temperature]', speed='[speed]')
def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
r"""Calculate the Wind Chill Temperature Index (WCTI).
r'''Calculate the Wind Chill Temperature Index (WCTI).

Calculates WCTI from the current temperature and wind speed using the formula
outlined by the FCM [FCMR192003]_.
Expand Down Expand Up @@ -154,7 +154,7 @@ def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
--------
heat_index

"""
'''
# Correct for lower height measurement of winds if necessary
if face_level_winds:
# No in-place so that we copy
Expand All @@ -178,7 +178,7 @@ def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
@exporter.export
@check_units('[temperature]')
def heat_index(temperature, rh, mask_undefined=True):
r"""Calculate the Heat Index from the current temperature and relative humidity.
r'''Calculate the Heat Index from the current temperature and relative humidity.

The implementation uses the formula outlined in [Rothfusz1990]_. This equation is a
multi-variable least-squares regression of the values obtained in [Steadman1979]_.
Expand Down Expand Up @@ -208,7 +208,7 @@ def heat_index(temperature, rh, mask_undefined=True):
--------
windchill

"""
'''
delta = temperature - 0. * units.degF
rh2 = rh * rh
delta2 = delta * delta
Expand All @@ -232,7 +232,7 @@ def heat_index(temperature, rh, mask_undefined=True):
@exporter.export
@check_units('[pressure]')
def pressure_to_height_std(pressure):
r"""Convert pressure data to heights using the U.S. standard atmosphere.
r'''Convert pressure data to heights using the U.S. standard atmosphere.

The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61.

Expand All @@ -250,17 +250,98 @@ def pressure_to_height_std(pressure):
-----
.. math:: Z = \frac{T_0}{\Gamma}[1-\frac{p}{p_0}^\frac{R\Gamma}{g}]

"""
'''
t0 = 288. * units.kelvin
gamma = 6.5 * units('K/km')
p0 = 1013.25 * units.mbar
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.

'''
# 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.

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

return height


@exporter.export
@check_units('[length]')
def height_to_pressure_std(height):
r"""Convert height data to pressures using the U.S. standard atmosphere.
r'''Convert height data to pressures using the U.S. standard atmosphere.

The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61.

Expand All @@ -278,7 +359,7 @@ def height_to_pressure_std(height):
-----
.. math:: p = p_0 e^{\frac{g}{R \Gamma} \text{ln}(1-\frac{Z \Gamma}{T_0})}

"""
'''
t0 = 288. * units.kelvin
gamma = 6.5 * units('K/km')
p0 = 1013.25 * units.mbar
Expand All @@ -287,7 +368,7 @@ def height_to_pressure_std(height):

@exporter.export
def coriolis_parameter(latitude):
r"""Calculate the coriolis parameter at each point.
r'''Calculate the coriolis parameter at each point.

The implementation uses the formula outlined in [Hobbs1977]_ pg.370-371.

Expand All @@ -301,15 +382,15 @@ def coriolis_parameter(latitude):
`pint.Quantity`
The corresponding coriolis force at each point

"""
'''
latitude = _check_radians(latitude, max_radians=np.pi / 2)
return 2. * omega * np.sin(latitude)


@exporter.export
@check_units('[pressure]', '[length]')
def add_height_to_pressure(pressure, height):
r"""Calculate the pressure at a certain height above another pressure level.
r'''Calculate the pressure at a certain height above another pressure level.

This assumes a standard atmosphere.

Expand All @@ -329,15 +410,15 @@ def add_height_to_pressure(pressure, height):
--------
pressure_to_height_std, height_to_pressure_std, add_pressure_to_height

"""
'''
pressure_level_height = pressure_to_height_std(pressure)
return height_to_pressure_std(pressure_level_height + height)


@exporter.export
@check_units('[length]', '[pressure]')
def add_pressure_to_height(height, pressure):
r"""Calculate the height at a certain pressure above another height.
r'''Calculate the height at a certain pressure above another height.

This assumes a standard atmosphere.

Expand All @@ -357,15 +438,15 @@ def add_pressure_to_height(height, pressure):
--------
pressure_to_height_std, height_to_pressure_std, add_height_to_pressure

"""
'''
pressure_at_height = height_to_pressure_std(height)
return pressure_to_height_std(pressure_at_height - pressure)


@exporter.export
@check_units('[dimensionless]', '[pressure]', '[pressure]')
def sigma_to_pressure(sigma, psfc, ptop):
r"""Calculate pressure from sigma values.
r'''Calculate pressure from sigma values.

Parameters
----------
Expand Down Expand Up @@ -394,7 +475,7 @@ def sigma_to_pressure(sigma, psfc, ptop):
* :math:`p_{sfc}` is pressure at the surface or model floor
* :math:`p_{top}` is pressure at the top of the model domain

"""
'''
if np.any(sigma < 0) or np.any(sigma > 1):
raise ValueError('Sigma values should be bounded by 0 and 1')

Expand All @@ -405,7 +486,7 @@ def sigma_to_pressure(sigma, psfc, ptop):


def _check_radians(value, max_radians=2 * np.pi):
"""Input validation of values that could be in degrees instead of radians.
'''Input validation of values that could be in degrees instead of radians.

Parameters
----------
Expand All @@ -420,7 +501,7 @@ def _check_radians(value, max_radians=2 * np.pi):
`pint.Quantity`
The input value

"""
'''
try:
value = value.to('radians').m
except AttributeError:
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 = gravit_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.

I think you probably meant "gravitational_constant" here.

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