# Copyright (c) 2008,2015,2016,2017,2018,2019 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains a collection of basic calculations.
These include:
* wind components
* heat index
* windchill
"""
import contextlib
from itertools import product
import warnings
import numpy as np
from scipy.ndimage import gaussian_filter
from .. import constants as mpconsts
from ..package_tools import Exporter
from ..units import check_units, masked_array, units
from ..xarray import preprocess_and_wrap
exporter = Exporter(globals())
# The following variables are constants for a standard atmosphere
t0 = units.Quantity(288., 'kelvin')
p0 = units.Quantity(1013.25, 'hPa')
gamma = units.Quantity(6.5, 'K/km')
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='u')
@check_units('[speed]', '[speed]')
def wind_speed(u, v):
r"""Compute the wind speed from u and v-components.
Parameters
----------
u : `pint.Quantity`
Wind component in the X (East-West) direction
v : `pint.Quantity`
Wind component in the Y (North-South) direction
Returns
-------
wind speed: `pint.Quantity`
Speed of the wind
See Also
--------
wind_components
"""
speed = np.sqrt(u * u + v * v)
return speed
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='u')
@check_units('[speed]', '[speed]')
def wind_direction(u, v, convention='from'):
r"""Compute the wind direction from u and v-components.
Parameters
----------
u : `pint.Quantity`
Wind component in the X (East-West) direction
v : `pint.Quantity`
Wind component in the Y (North-South) direction
convention : str
Convention to return direction; 'from' returns the direction the wind is coming from
(meteorological convention), 'to' returns the direction the wind is going towards
(oceanographic convention), default is 'from'.
Returns
-------
direction: `pint.Quantity`
The direction of the wind in intervals [0, 360] degrees, with 360 being North,
direction defined by the convention kwarg.
See Also
--------
wind_components
Notes
-----
In the case of calm winds (where `u` and `v` are zero), this function returns a direction
of 0.
"""
wdir = units.Quantity(90., 'deg') - np.arctan2(-v, -u)
origshape = wdir.shape
wdir = np.atleast_1d(wdir)
# Handle oceanographic convection
if convention == 'to':
wdir -= units.Quantity(180., 'deg')
elif convention not in ('to', 'from'):
raise ValueError('Invalid kwarg for "convention". Valid options are "from" or "to".')
mask = wdir <= 0
if np.any(mask):
wdir[mask] += units.Quantity(360., 'deg')
# avoid unintended modification of `pint.Quantity` by direct use of magnitude
calm_mask = (np.asarray(u.magnitude) == 0.) & (np.asarray(v.magnitude) == 0.)
# np.any check required for legacy numpy which treats 0-d False boolean index as zero
if np.any(calm_mask):
wdir[calm_mask] = units.Quantity(0., 'deg')
return wdir.reshape(origshape).to('degrees')
[docs]@exporter.export
@preprocess_and_wrap(wrap_like=('speed', 'speed'))
@check_units('[speed]')
def wind_components(speed, wind_direction):
r"""Calculate the U, V wind vector components from the speed and direction.
Parameters
----------
speed : `pint.Quantity`
Wind speed (magnitude)
wind_direction : `pint.Quantity`
Wind direction, specified as the direction from which the wind is
blowing (0-2 pi radians or 0-360 degrees), with 360 degrees being North.
Returns
-------
u, v : tuple of `pint.Quantity`
The wind components in the X (East-West) and Y (North-South)
directions, respectively.
See Also
--------
wind_speed
wind_direction
Examples
--------
>>> from metpy.units import units
>>> metpy.calc.wind_components(10. * units('m/s'), 225. * units.deg)
(<Quantity(7.07106781, 'meter / second')>, <Quantity(7.07106781, 'meter / second')>)
.. versionchanged:: 1.0
Renamed ``wdir`` parameter to ``wind_direction``
"""
wind_direction = _check_radians(wind_direction, max_radians=4 * np.pi)
u = -speed * np.sin(wind_direction)
v = -speed * np.cos(wind_direction)
return u, v
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='temperature')
@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).
Calculates WCTI from the current temperature and wind speed using the formula
outlined by the FCM [FCMR192003]_.
Specifically, these formulas assume that wind speed is measured at
10m. If, instead, the speeds are measured at face level, the winds
need to be multiplied by a factor of 1.5 (this can be done by specifying
`face_level_winds` as `True`).
Parameters
----------
temperature : `pint.Quantity`
Air temperature
speed : `pint.Quantity`
Wind speed at 10m. If instead the winds are at face level,
`face_level_winds` should be set to `True` and the 1.5 multiplicative
correction will be applied automatically.
face_level_winds : bool, optional
A flag indicating whether the wind speeds were measured at facial
level instead of 10m, thus requiring a correction. Defaults to
`False`.
mask_undefined : bool, optional
A flag indicating whether a masked array should be returned with
values where wind chill is undefined masked. These are values where
the temperature > 50F or wind speed <= 3 miles per hour. Defaults
to `True`.
Returns
-------
`pint.Quantity`
Corresponding Wind Chill Temperature Index value(s)
See Also
--------
heat_index, apparent_temperature
"""
# Correct for lower height measurement of winds if necessary
if face_level_winds:
# No in-place so that we copy
# noinspection PyAugmentAssignment
speed = speed * 1.5
temp_limit, speed_limit = units.Quantity(10., 'degC'), units.Quantity(3, 'mph')
speed_factor = speed.to('km/hr').magnitude ** 0.16
wcti = units.Quantity((0.6215 + 0.3965 * speed_factor) * temperature.to('degC').magnitude
- 11.37 * speed_factor + 13.12, units.degC).to(temperature.units)
# See if we need to mask any undefined values
if mask_undefined:
mask = np.array((temperature > temp_limit) | (speed <= speed_limit))
if mask.any():
wcti = masked_array(wcti, mask=mask)
return wcti
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='temperature')
@check_units('[temperature]')
def heat_index(temperature, relative_humidity, mask_undefined=True):
r"""Calculate the Heat Index from the current temperature and relative humidity.
The implementation uses the formula outlined in [Rothfusz1990]_, which is a
multi-variable least-squares regression of the values obtained in [Steadman1979]_.
Additional conditional corrections are applied to match what the National
Weather Service operationally uses. See Figure 3 of [Anderson2013]_ for a
depiction of this algorithm and further discussion.
Parameters
----------
temperature : `pint.Quantity`
Air temperature
relative_humidity : `pint.Quantity`
The relative humidity expressed as a unitless ratio in the range [0, 1].
Can also pass a percentage if proper units are attached.
Returns
-------
`pint.Quantity`
Corresponding Heat Index value(s)
Other Parameters
----------------
mask_undefined : bool, optional
A flag indicating whether a masked array should be returned with
values masked where the temperature < 80F. Defaults to `True`.
.. versionchanged:: 1.0
Renamed ``rh`` parameter to ``relative_humidity``
See Also
--------
windchill, apparent_temperature
"""
temperature = np.atleast_1d(temperature)
relative_humidity = np.atleast_1d(relative_humidity)
# assign units to relative_humidity if they currently are not present
if not hasattr(relative_humidity, 'units'):
relative_humidity = units.Quantity(relative_humidity, 'dimensionless')
delta = temperature.to(units.degF) - units.Quantity(0., 'degF')
rh2 = relative_humidity * relative_humidity
delta2 = delta * delta
# Simplifed Heat Index -- constants converted for relative_humidity in [0, 1]
a = (units.Quantity(-10.3, 'degF') + 1.1 * delta
+ units.Quantity(4.7, 'delta_degF') * relative_humidity)
# More refined Heat Index -- constants converted for relative_humidity in [0, 1]
b = (units.Quantity(-42.379, 'degF')
+ 2.04901523 * delta
+ units.Quantity(1014.333127, 'delta_degF') * relative_humidity
- 22.475541 * delta * relative_humidity
- units.Quantity(6.83783e-3, '1/delta_degF') * delta2
- units.Quantity(5.481717e2, 'delta_degF') * rh2
+ units.Quantity(1.22874e-1, '1/delta_degF') * delta2 * relative_humidity
+ 8.5282 * delta * rh2
- units.Quantity(1.99e-2, '1/delta_degF') * delta2 * rh2)
# Create return heat index
hi = units.Quantity(np.full(np.shape(temperature), np.nan), 'degF')
# Retain masked status of temperature with resulting heat index
if hasattr(temperature, 'mask'):
hi = masked_array(hi)
# If T <= 40F, Heat Index is T
sel = (temperature <= units.Quantity(40., 'degF'))
if np.any(sel):
hi[sel] = temperature[sel].to(units.degF)
# If a < 79F and hi is unset, Heat Index is a
sel = (a < units.Quantity(79., 'degF')) & np.isnan(hi)
if np.any(sel):
hi[sel] = a[sel]
# Use b now for anywhere hi has yet to be set
sel = np.isnan(hi)
if np.any(sel):
hi[sel] = b[sel]
# Adjustment for RH <= 13% and 80F <= T <= 112F
sel = ((relative_humidity <= units.Quantity(13., 'percent'))
& (temperature >= units.Quantity(80., 'degF'))
& (temperature <= units.Quantity(112., 'degF')))
if np.any(sel):
rh15adj = ((13. - relative_humidity[sel] * 100.) / 4.
* np.sqrt((units.Quantity(17., 'delta_degF')
- np.abs(delta[sel] - units.Quantity(95., 'delta_degF')))
/ units.Quantity(17., 'delta_degF')))
hi[sel] = hi[sel] - units.Quantity(rh15adj, 'delta_degF')
# Adjustment for RH > 85% and 80F <= T <= 87F
sel = ((relative_humidity > units.Quantity(85., 'percent'))
& (temperature >= units.Quantity(80., 'degF'))
& (temperature <= units.Quantity(87., 'degF')))
if np.any(sel):
rh85adj = (0.02 * (relative_humidity[sel] * 100. - 85.)
* (units.Quantity(87., 'delta_degF') - delta[sel]))
hi[sel] = hi[sel] + rh85adj
# See if we need to mask any undefined values
if mask_undefined:
mask = np.array(temperature < units.Quantity(80., 'degF'))
if mask.any():
hi = masked_array(hi, mask=mask)
return hi
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='temperature')
@check_units(temperature='[temperature]', speed='[speed]')
def apparent_temperature(temperature, relative_humidity, speed, face_level_winds=False,
mask_undefined=True):
r"""Calculate the current apparent temperature.
Calculates the current apparent temperature based on the wind chill or heat index
as appropriate for the current conditions. Follows [NWS10201]_.
Parameters
----------
temperature : `pint.Quantity`
Air temperature
relative_humidity : `pint.Quantity`
Relative humidity expressed as a unitless ratio in the range [0, 1].
Can also pass a percentage if proper units are attached.
speed : `pint.Quantity`
Wind speed at 10m. If instead the winds are at face level,
`face_level_winds` should be set to `True` and the 1.5 multiplicative
correction will be applied automatically.
face_level_winds : bool, optional
A flag indicating whether the wind speeds were measured at facial
level instead of 10m, thus requiring a correction. Defaults to
`False`.
mask_undefined : bool, optional
A flag indicating whether a masked array should be returned with
values where wind chill or heat_index is undefined masked. For wind
chill, these are values where the temperature > 50F or
wind speed <= 3 miles per hour. For heat index, these are values
where the temperature < 80F.
Defaults to `True`.
Returns
-------
`pint.Quantity`
Corresponding apparent temperature value(s)
.. versionchanged:: 1.0
Renamed ``rh`` parameter to ``relative_humidity``
See Also
--------
heat_index, windchill
"""
is_not_scalar = isinstance(temperature.m, (list, tuple, np.ndarray))
temperature = np.atleast_1d(temperature)
relative_humidity = np.atleast_1d(relative_humidity)
speed = np.atleast_1d(speed)
# NB: mask_defined=True is needed to know where computed values exist
wind_chill_temperature = windchill(temperature, speed, face_level_winds=face_level_winds,
mask_undefined=True).to(temperature.units)
heat_index_temperature = heat_index(temperature, relative_humidity,
mask_undefined=True).to(temperature.units)
# Combine the heat index and wind chill arrays (no point has a value in both)
# NB: older numpy.ma.where does not return a masked array
app_temperature = masked_array(
np.ma.where(masked_array(wind_chill_temperature).mask,
heat_index_temperature.to(temperature.units),
wind_chill_temperature.to(temperature.units)
), temperature.units)
# If mask_undefined is False, then set any masked values to the temperature
if not mask_undefined:
app_temperature[app_temperature.mask] = temperature[app_temperature.mask]
# If no values are masked and provided temperature does not have a mask
# we should return a non-masked array
if not np.any(app_temperature.mask) and not hasattr(temperature, 'mask'):
app_temperature = units.Quantity(np.array(app_temperature.m), temperature.units)
if is_not_scalar:
return app_temperature
else:
return np.atleast_1d(app_temperature)[0]
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='pressure')
@check_units('[pressure]')
def pressure_to_height_std(pressure):
r"""Convert pressure data to height using the U.S. standard atmosphere [NOAA1976]_.
The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61.
Parameters
----------
pressure : `pint.Quantity`
Atmospheric pressure
Returns
-------
`pint.Quantity`
Corresponding height value(s)
Notes
-----
.. math:: Z = \frac{T_0}{\Gamma}[1-\frac{p}{p_0}^\frac{R\Gamma}{g}]
"""
return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(
mpconsts.Rd * gamma / mpconsts.g))
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='height')
@check_units('[length]')
def height_to_geopotential(height):
r"""Compute geopotential for a given height above sea level.
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 = \frac{g R_e z}{R_e + z}
(where :math:`\Phi` is geopotential, :math:`z` is height, :math:`R_e` is average Earth
radius, and :math:`g` is standard gravity).
Parameters
----------
height : `pint.Quantity`
Height above sea level
Returns
-------
`pint.Quantity`
Corresponding geopotential value(s)
Examples
--------
>>> 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. 9805.11097983 19607.1448853 29406.10316465
39201.98726524 48994.79863351 58784.53871501 68571.20895435
78354.81079527 88135.34568058 97912.81505219], '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
"""
return (mpconsts.g * mpconsts.Re * height) / (mpconsts.Re + height)
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='geopotential')
def geopotential_to_height(geopotential):
r"""Compute height above sea level from a given geopotential.
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{\Phi R_e}{gR_e - \Phi}
(where :math:`\Phi` is geopotential, :math:`z` is height, :math:`R_e` is average Earth
radius, and :math:`g` is standard gravity).
Parameters
----------
geopotential : `pint.Quantity`
Geopotential
Returns
-------
`pint.Quantity`
Corresponding value(s) of height above sea level
Examples
--------
>>> 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. 9805.11097983 19607.1448853 29406.10316465
39201.98726524 48994.79863351 58784.53871501 68571.20895435
78354.81079527 88135.34568058 97912.81505219], '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.)
.. versionchanged:: 1.0
Renamed ``geopot`` parameter to ``geopotential``
See Also
--------
height_to_geopotential
"""
return (geopotential * mpconsts.Re) / (mpconsts.g * mpconsts.Re - geopotential)
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='height')
@check_units('[length]')
def height_to_pressure_std(height):
r"""Convert height data to pressures using the U.S. standard atmosphere [NOAA1976]_.
The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61.
Parameters
----------
height : `pint.Quantity`
Atmospheric height
Returns
-------
`pint.Quantity`
Corresponding pressure value(s)
Notes
-----
.. math:: p = p_0 e^{\frac{g}{R \Gamma} \text{ln}(1-\frac{Z \Gamma}{T_0})}
"""
return p0 * (1 - (gamma / t0) * height) ** (mpconsts.g / (mpconsts.Rd * gamma))
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='latitude')
def coriolis_parameter(latitude):
r"""Calculate the coriolis parameter at each point.
The implementation uses the formula outlined in [Hobbs1977]_ pg.370-371.
Parameters
----------
latitude : array_like
Latitude at each point
Returns
-------
`pint.Quantity`
Corresponding coriolis force at each point
"""
latitude = _check_radians(latitude, max_radians=np.pi / 2)
return (2. * mpconsts.omega * np.sin(latitude)).to('1/s')
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='pressure')
@check_units('[pressure]', '[length]')
def add_height_to_pressure(pressure, height):
r"""Calculate the pressure at a certain height above another pressure level.
This assumes a standard atmosphere [NOAA1976]_.
Parameters
----------
pressure : `pint.Quantity`
Pressure level
height : `pint.Quantity`
Height above a pressure level
Returns
-------
`pint.Quantity`
Corresponding pressure value for the height above the pressure level
See Also
--------
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)
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='height')
@check_units('[length]', '[pressure]')
def add_pressure_to_height(height, pressure):
r"""Calculate the height at a certain pressure above another height.
This assumes a standard atmosphere [NOAA1976]_.
Parameters
----------
height : `pint.Quantity`
Height level
pressure : `pint.Quantity`
Pressure above height level
Returns
-------
`pint.Quantity`
The corresponding height value for the pressure above the height level
See Also
--------
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)
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='sigma')
@check_units('[dimensionless]', '[pressure]', '[pressure]')
def sigma_to_pressure(sigma, pressure_sfc, pressure_top):
r"""Calculate pressure from sigma values.
Parameters
----------
sigma : ndarray
Sigma levels to be converted to pressure levels
pressure_sfc : `pint.Quantity`
Surface pressure value
pressure_top : `pint.Quantity`
Pressure value at the top of the model domain
Returns
-------
`pint.Quantity`
Pressure values at the given sigma levels
Notes
-----
Sigma definition adapted from [Philips1957]_:
.. math:: p = \sigma * (p_{sfc} - p_{top}) + p_{top}
* :math:`p` is pressure at a given `\sigma` level
* :math:`\sigma` is non-dimensional, scaled pressure
* :math:`p_{sfc}` is pressure at the surface or model floor
* :math:`p_{top}` is pressure at the top of the model domain
.. versionchanged:: 1.0
Renamed ``psfc``, ``ptop`` parameters to ``pressure_sfc``, ``pressure_top``
"""
if np.any(sigma < 0) or np.any(sigma > 1):
raise ValueError('Sigma values should be bounded by 0 and 1')
if pressure_sfc.magnitude < 0 or pressure_top.magnitude < 0:
raise ValueError('Pressure input should be non-negative')
return sigma * (pressure_sfc - pressure_top) + pressure_top
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='scalar_grid', match_unit=True, to_magnitude=True)
def smooth_gaussian(scalar_grid, n):
"""Filter with normal distribution of weights.
Parameters
----------
scalar_grid : `pint.Quantity`
Some n-dimensional scalar grid. If more than two axes, smoothing
is only done across the last two.
n : int
Degree of filtering
Returns
-------
`pint.Quantity`
The filtered 2D scalar grid
Notes
-----
This function is a close replication of the GEMPAK function ``GWFS``,
but is not identical. The following notes are incorporated from
the GEMPAK source code:
This function smoothes a scalar grid using a moving average
low-pass filter whose weights are determined by the normal
(Gaussian) probability distribution function for two dimensions.
The weight given to any grid point within the area covered by the
moving average for a target grid point is proportional to:
.. math:: e^{-D^2}
where D is the distance from that point to the target point divided
by the standard deviation of the normal distribution. The value of
the standard deviation is determined by the degree of filtering
requested. The degree of filtering is specified by an integer.
This integer is the number of grid increments from crest to crest
of the wave for which the theoretical response is 1/e = .3679. If
the grid increment is called delta_x, and the value of this integer
is represented by N, then the theoretical filter response function
value for the N * delta_x wave will be 1/e. The actual response
function will be greater than the theoretical value.
The larger N is, the more severe the filtering will be, because the
response function for all wavelengths shorter than N * delta_x
will be less than 1/e. Furthermore, as N is increased, the slope
of the filter response function becomes more shallow; so, the
response at all wavelengths decreases, but the amount of decrease
lessens with increasing wavelength. (The theoretical response
function can be obtained easily--it is the Fourier transform of the
weight function described above.)
The area of the patch covered by the moving average varies with N.
As N gets bigger, the smoothing gets stronger, and weight values
farther from the target grid point are larger because the standard
deviation of the normal distribution is bigger. Thus, increasing
N has the effect of expanding the moving average window as well as
changing the values of weights. The patch is a square covering all
points whose weight values are within two standard deviations of the
mean of the two dimensional normal distribution.
The key difference between GEMPAK's GWFS and this function is that,
in GEMPAK, the leftover weight values representing the fringe of the
distribution are applied to the target grid point. In this
function, the leftover weights are not used.
When this function is invoked, the first argument is the grid to be
smoothed, the second is the value of N as described above:
GWFS ( S, N )
where N > 1. If N <= 1, N = 2 is assumed. For example, if N = 4,
then the 4 delta x wave length is passed with approximate response
1/e.
"""
# Compute standard deviation in a manner consistent with GEMPAK
n = int(round(n))
if n < 2:
n = 2
sgma = n / (2 * np.pi)
# Construct sigma sequence so smoothing occurs only in horizontal direction
nax = len(scalar_grid.shape)
# Assume the last two axes represent the horizontal directions
sgma_seq = [sgma if i > nax - 3 else 0 for i in range(nax)]
# Compute smoothed field
return gaussian_filter(scalar_grid, sgma_seq, truncate=2 * np.sqrt(2))
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='scalar_grid', match_unit=True, to_magnitude=True)
def smooth_window(scalar_grid, window, passes=1, normalize_weights=True):
"""Filter with an arbitrary window smoother.
Parameters
----------
scalar_grid : array-like
N-dimensional scalar grid to be smoothed
window : ndarray
Window to use in smoothing. Can have dimension less than or equal to N. If
dimension less than N, the scalar grid will be smoothed along its trailing dimensions.
Shape along each dimension must be odd.
passes : int
The number of times to apply the filter to the grid. Defaults to 1.
normalize_weights : bool
If true, divide the values in window by the sum of all values in the window to obtain
the normalized smoothing weights. If false, use supplied values directly as the
weights.
Returns
-------
array-like
The filtered scalar grid
Notes
-----
This function can be applied multiple times to create a more smoothed field and will only
smooth the interior points, leaving the end points with their original values (this
function will leave an unsmoothed edge of size `(n - 1) / 2` for each `n` in the shape of
`window` around the data). If a masked value or NaN values exists in the array, it will
propagate to any point that uses that particular grid point in the smoothing calculation.
Applying the smoothing function multiple times will propogate NaNs further throughout the
domain.
See Also
--------
smooth_rectangular, smooth_circular, smooth_n_point, smooth_gaussian
"""
def _pad(n):
# Return number of entries to pad given length along dimension.
return (n - 1) // 2
def _zero_to_none(x):
# Convert zero values to None, otherwise return what is given.
return x if x != 0 else None
def _offset(pad, k):
# Return padded slice offset by k entries
return slice(_zero_to_none(pad + k), _zero_to_none(-pad + k))
def _trailing_dims(indexer):
# Add ... to the front of an indexer, since we are working with trailing dimensions.
return (Ellipsis,) + tuple(indexer)
# Verify that shape in all dimensions is odd (need to have a neighboorhood around a
# central point)
if any((size % 2 == 0) for size in window.shape):
raise ValueError('The shape of the smoothing window must be odd in all dimensions.')
# Optionally normalize the supplied weighting window
if normalize_weights:
weights = window / np.sum(window)
else:
weights = window
# Set indexes
# Inner index for the centered array elements that are affected by the smoothing
inner_full_index = _trailing_dims(_offset(_pad(n), 0) for n in weights.shape)
# Indexes to iterate over each weight
weight_indexes = tuple(product(*(range(n) for n in weights.shape)))
# Index for full array elements, offset by the weight index
def offset_full_index(weight_index):
return _trailing_dims(_offset(_pad(n), weight_index[i] - _pad(n))
for i, n in enumerate(weights.shape))
# TODO: this is not lazy-loading/dask compatible, as it "densifies" the data
data = np.array(scalar_grid)
for _ in range(passes):
# Set values corresponding to smoothing weights by summing over each weight and
# applying offsets in needed dimensions
data[inner_full_index] = sum(weights[index] * data[offset_full_index(index)]
for index in weight_indexes)
return data
[docs]@exporter.export
def smooth_rectangular(scalar_grid, size, passes=1):
"""Filter with a rectangular window smoother.
Parameters
----------
scalar_grid : array-like
N-dimensional scalar grid to be smoothed
size : int or sequence of ints
Shape of rectangle along the trailing dimension(s) of the scalar grid
passes : int
The number of times to apply the filter to the grid. Defaults to 1.
Returns
-------
array-like
The filtered scalar grid
Notes
-----
This function can be applied multiple times to create a more smoothed field and will only
smooth the interior points, leaving the end points with their original values (this
function will leave an unsmoothed edge of size `(n - 1) / 2` for each `n` in `size` around
the data). If a masked value or NaN values exists in the array, it will propagate to any
point that uses that particular grid point in the smoothing calculation. Applying the
smoothing function multiple times will propogate NaNs further throughout the domain.
See Also
--------
smooth_window, smooth_circular, smooth_n_point, smooth_gaussian
"""
return smooth_window(scalar_grid, np.ones(size), passes=passes)
[docs]@exporter.export
def smooth_circular(scalar_grid, radius, passes=1):
"""Filter with a circular window smoother.
Parameters
----------
scalar_grid : array-like
N-dimensional scalar grid to be smoothed. If more than two axes, smoothing is only
done along the last two.
radius : int
Radius of the circular smoothing window. The "diameter" of the circle (width of
smoothing window) is 2 * radius + 1 to provide a smoothing window with odd shape.
passes : int
The number of times to apply the filter to the grid. Defaults to 1.
Returns
-------
array-like
The filtered scalar grid
Notes
-----
This function can be applied multiple times to create a more smoothed field and will only
smooth the interior points, leaving the end points with their original values (this
function will leave an unsmoothed edge of size `radius` around the data). If a masked
value or NaN values exists in the array, it will propagate to any point that uses that
particular grid point in the smoothing calculation. Applying the smoothing function
multiple times will propogate NaNs further throughout the domain.
See Also
--------
smooth_window, smooth_rectangular, smooth_n_point, smooth_gaussian
"""
# Generate the circle
size = 2 * radius + 1
x, y = np.mgrid[:size, :size]
distance = np.sqrt((x - radius) ** 2 + (y - radius) ** 2)
circle = distance <= radius
# Apply smoother
return smooth_window(scalar_grid, circle, passes=passes)
[docs]@exporter.export
def smooth_n_point(scalar_grid, n=5, passes=1):
"""Filter with an n-point smoother.
Parameters
----------
scalar_grid : array-like or `pint.Quantity`
N-dimensional scalar grid to be smoothed. If more than two axes, smoothing is only
done along the last two.
n: int
The number of points to use in smoothing, only valid inputs
are 5 and 9. Defaults to 5.
passes : int
The number of times to apply the filter to the grid. Defaults to 1.
Returns
-------
array-like or `pint.Quantity`
The filtered scalar grid
Notes
-----
This function is a close replication of the GEMPAK function SM5S and SM9S depending on the
choice of the number of points to use for smoothing. This function can be applied multiple
times to create a more smoothed field and will only smooth the interior points, leaving
the end points with their original values (this function will leave an unsmoothed edge of
size 1 around the data). If a masked value or NaN values exists in the array, it will
propagate to any point that uses that particular grid point in the smoothing calculation.
Applying the smoothing function multiple times will propogate NaNs further throughout the
domain.
See Also
--------
smooth_window, smooth_rectangular, smooth_circular, smooth_gaussian
"""
if n == 9:
weights = np.array([[0.0625, 0.125, 0.0625],
[0.125, 0.25, 0.125],
[0.0625, 0.125, 0.0625]])
elif n == 5:
weights = np.array([[0., 0.125, 0.],
[0.125, 0.5, 0.125],
[0., 0.125, 0.]])
else:
raise ValueError('The number of points to use in the smoothing '
'calculation must be either 5 or 9.')
return smooth_window(scalar_grid, window=weights, passes=passes, normalize_weights=False)
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='altimeter_value')
@check_units('[pressure]', '[length]')
def altimeter_to_station_pressure(altimeter_value, height):
r"""Convert the altimeter measurement to station pressure.
This function is useful for working with METARs since they do not provide
altimeter values, but not sea-level pressure or station pressure.
The following definitions of altimeter setting and station pressure
are taken from [Smithsonian1951]_ Altimeter setting is the
pressure value to which an aircraft altimeter scale is set so that it will
indicate the altitude above mean sea-level of an aircraft on the ground at the
location for which the value is determined. It assumes a standard atmosphere [NOAA1976]_.
Station pressure is the atmospheric pressure at the designated station elevation.
Finding the station pressure can be helpful for calculating sea-level pressure
or other parameters.
Parameters
----------
altimeter_value : `pint.Quantity`
The altimeter setting value as defined by the METAR or other observation,
which can be measured in either inches of mercury (in. Hg) or millibars (mb)
height: `pint.Quantity`
Elevation of the station measuring pressure
Returns
-------
`pint.Quantity`
The station pressure in hPa or in. Hg. Can be used to calculate sea-level
pressure.
See Also
--------
altimeter_to_sea_level_pressure
Notes
-----
This function is implemented using the following equations from the
Smithsonian Handbook (1951) p. 269
Equation 1:
.. math:: A_{mb} = (p_{mb} - 0.3)F
Equation 3:
.. math:: F = \left [1 + \left(\frac{p_{0}^n a}{T_{0}} \right)
\frac{H_{b}}{p_{1}^n} \right ] ^ \frac{1}{n}
Where,
:math:`p_{0}` = standard sea-level pressure = 1013.25 mb
:math:`p_{1} = p_{mb} - 0.3` when :math:`p_{0} = 1013.25 mb`
gamma = lapse rate in [NOAA1976]_ standard atmosphere below the isothermal layer
:math:`6.5^{\circ}C. km.^{-1}`
:math:`t_{0}` = standard sea-level temperature 288 K
:math:`H_{b} =` station elevation in meters (elevation for which station pressure is given)
:math:`n = \frac{a R_{d}}{g} = 0.190284` where :math:`R_{d}` is the gas constant for dry
air
And solving for :math:`p_{mb}` results in the equation below, which is used to
calculate station pressure :math:`(p_{mb})`
.. math:: p_{mb} = \left [A_{mb} ^ n - \left (\frac{p_{0} a H_{b}}{T_0}
\right) \right] ^ \frac{1}{n} + 0.3
"""
# N-Value
n = (mpconsts.Rd * gamma / mpconsts.g).to_base_units()
return ((altimeter_value ** n
- ((p0.to(altimeter_value.units) ** n * gamma * height) / t0)) ** (1 / n)
+ units.Quantity(0.3, 'hPa'))
[docs]@exporter.export
@preprocess_and_wrap(wrap_like='altimeter_value')
@check_units('[pressure]', '[length]', '[temperature]')
def altimeter_to_sea_level_pressure(altimeter_value, height, temperature):
r"""Convert the altimeter setting to sea-level pressure.
This function is useful for working with METARs since most provide
altimeter values, but not sea-level pressure, which is often plotted
on surface maps. The following definitions of altimeter setting, station pressure, and
sea-level pressure are taken from [Smithsonian1951]_.
Altimeter setting is the pressure value to which an aircraft altimeter scale
is set so that it will indicate the altitude above mean sea-level of an aircraft
on the ground at the location for which the value is determined. It assumes a standard
atmosphere. Station pressure is the atmospheric pressure at the designated station
elevation. Sea-level pressure is a pressure value obtained by the theoretical reduction
of barometric pressure to sea level. It is assumed that atmosphere extends to sea level
below the station and that the properties of the atmosphere are related to conditions
observed at the station. This value is recorded by some surface observation stations,
but not all. If the value is recorded, it can be found in the remarks section. Finding
the sea-level pressure is helpful for plotting purposes and different calculations.
Parameters
----------
altimeter_value : 'pint.Quantity'
The altimeter setting value is defined by the METAR or other observation,
with units of inches of mercury (in Hg) or millibars (hPa).
height : 'pint.Quantity'
Elevation of the station measuring pressure. Often times measured in meters
temperature : 'pint.Quantity'
Temperature at the station
Returns
-------
'pint.Quantity'
The sea-level pressure in hPa and makes pressure values easier to compare
between different stations.
See Also
--------
altimeter_to_station_pressure
Notes
-----
This function is implemented using the following equations from Wallace and Hobbs (1977).
Equation 2.29:
.. math::
\Delta z = Z_{2} - Z_{1}
= \frac{R_{d} \bar T_{v}}{g_0}ln\left(\frac{p_{1}}{p_{2}}\right)
= \bar H ln \left (\frac {p_{1}}{p_{2}} \right)
Equation 2.31:
.. math::
p_{0} = p_{g}exp \left(\frac{Z_{g}}{\bar H} \right) \\
= p_{g}exp \left(\frac{g_{0}Z_{g}}{R_{d}\bar T_{v}} \right)
Then by substituting :math:`Delta_{Z}` for :math:`Z_{g}` in Equation 2.31:
.. math:: p_{sea_level} = p_{station} exp\left(\frac{\Delta z}{H}\right)
where :math:`Delta_{Z}` is the elevation in meters and :math:`H = \frac{R_{d}T}{g}`
"""
# Calculate the station pressure using function altimeter_to_station_pressure()
psfc = altimeter_to_station_pressure(altimeter_value, height)
# Calculate the scale height
h = mpconsts.Rd * temperature / mpconsts.g
return psfc * np.exp(height / h)
def _check_radians(value, max_radians=2 * np.pi):
"""Input validation of values that could be in degrees instead of radians.
Parameters
----------
value : `pint.Quantity`
Input value to check
max_radians : float
Maximum absolute value of radians before warning
Returns
-------
`pint.Quantity`
Input value
"""
with contextlib.suppress(AttributeError):
value = value.to('radians').m
if np.any(np.greater(np.abs(value), max_radians)):
warnings.warn('Input over {} radians. '
'Ensure proper units are given.'.format(np.nanmax(max_radians)))
return value