Skip to content

Commit

Permalink
MNT: Move functions from future to their proper location (Fixes #1112)
Browse files Browse the repository at this point in the history
This makes their behavior standard and eliminates the old behavior.
  • Loading branch information
dopplershift committed Jan 10, 2020
1 parent 706c9f4 commit 67eeec7
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 220 deletions.
42 changes: 29 additions & 13 deletions src/metpy/calc/indices.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,11 @@
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains calculation of various derived indices."""
import warnings

import numpy as np

from .tools import get_layer
from .thermo import mixing_ratio, saturation_vapor_pressure
from .tools import _remove_nans, get_layer
from .. import constants as mpconsts
from ..package_tools import Exporter
from ..units import atleast_1d, check_units, concatenate, units
from ..xarray import preprocess_xarray
Expand All @@ -16,8 +16,8 @@

@exporter.export
@preprocess_xarray
@check_units('[temperature]', '[pressure]', '[pressure]')
def precipitable_water(dewpt, pressure, bottom=None, top=None):
@check_units('[pressure]', '[temperature]', bottom='[pressure]', top='[pressure]')
def precipitable_water(pressure, dewpt, *, bottom=None, top=None):
r"""Calculate precipitable water through the depth of a sounding.
Formula used is:
Expand All @@ -28,10 +28,10 @@ def precipitable_water(dewpt, pressure, bottom=None, top=None):
Parameters
----------
dewpt : `pint.Quantity`
Atmospheric dewpoint profile
pressure : `pint.Quantity`
Atmospheric pressure profile
dewpt : `pint.Quantity`
Atmospheric dewpoint profile
bottom: `pint.Quantity`, optional
Bottom of the layer, specified in pressure. Defaults to None (highest pressure).
top: `pint.Quantity`, optional
Expand All @@ -46,14 +46,30 @@ def precipitable_water(dewpt, pressure, bottom=None, top=None):
--------
>>> pressure = np.array([1000, 950, 900]) * units.hPa
>>> dewpoint = np.array([20, 15, 10]) * units.degC
>>> pw = precipitable_water(dewpoint, pressure)
>>> pw = precipitable_water(pressure, dewpoint)
"""
warnings.warn('Input variables will be reordered in 1.0 to be (pressure, dewpt, bottom,'
'top). To update to new input format before 1.0 is released, use'
'`from metpy.future import precipitable_water`.', FutureWarning)
from ..future import precipitable_water as _precipitable_water
return _precipitable_water(pressure, dewpt, bottom=bottom, top=top)
# Sort pressure and dewpoint to be in decreasing pressure order (increasing height)
sort_inds = np.argsort(pressure)[::-1]
pressure = pressure[sort_inds]
dewpt = dewpt[sort_inds]

pressure, dewpt = _remove_nans(pressure, dewpt)

if top is None:
top = np.nanmin(pressure.magnitude) * pressure.units

if bottom is None:
bottom = np.nanmax(pressure.magnitude) * pressure.units

pres_layer, dewpt_layer = get_layer(pressure, dewpt, bottom=bottom, depth=bottom - top)

w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)

# Since pressure is in decreasing order, pw will be the opposite sign of that expected.
pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units)
/ (mpconsts.g * mpconsts.rho_l))
return pw.to('millimeters')


@exporter.export
Expand Down
64 changes: 38 additions & 26 deletions src/metpy/calc/kinematics.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,11 @@
# SPDX-License-Identifier: BSD-3-Clause
"""Contains calculation of kinematic parameters (e.g. divergence or vorticity)."""
import functools
import warnings

import numpy as np

from . import coriolis_parameter
from .tools import first_derivative, gradient
from .tools import first_derivative, get_layer_heights, gradient
from .. import constants as mpconsts
from ..cbook import iterable
from ..package_tools import Exporter
Expand Down Expand Up @@ -447,14 +446,20 @@ def geostrophic_wind(heights, f, dx, dy):


@exporter.export
@check_units(f='[frequency]', dx='[length]', dy='[length]', u='[speed]', v='[speed]')
def ageostrophic_wind(heights, f, dx, dy, u, v, dim_order='yx'):
@preprocess_xarray
@ensure_yx_order
@check_units(f='[frequency]', u='[speed]', v='[speed]', dx='[length]', dy='[length]')
def ageostrophic_wind(heights, u, v, f, dx, dy, dim_order='yx'):
r"""Calculate the ageostrophic wind given from the heights or geopotential.
Parameters
----------
heights : (M, N) ndarray
The height or geopotential field.
u : (M, N) `pint.Quantity`
The u wind field.
v : (M, N) `pint.Quantity`
The u wind field.
f : array_like
The coriolis parameter. This can be a scalar to be applied
everywhere or an array of values.
Expand All @@ -464,30 +469,24 @@ def ageostrophic_wind(heights, f, dx, dy, u, v, dim_order='yx'):
dy : `pint.Quantity`
The grid spacing(s) in the y-direction. If an array, there should be one item less than
the size of `heights` along the applicable axis.
u : (M, N) `pint.Quantity`
The u wind field.
v : (M, N) `pint.Quantity`
The u wind field.
Returns
-------
A 2-item tuple of arrays, `pint.Quantity`
A 2-item tuple of arrays
A tuple of the u-component and v-component of the ageostrophic wind.
Notes
-----
If inputs have more than two dimensions, they are assumed to have either leading dimensions
of (x, y) or trailing dimensions of (y, x), depending on the value of ``dim_order``.
The order of the inputs will be changed in 1.0 to be (heights, u, v, f, dx, dy).
To updated to the new format, use `from metpy.future import ageostrophic_wind`.
This function contains an updated input variable order from the same function in the
kinematics module. This version will be fully implemented in 1.0 and moved from the
`future` module back to the `kinematics` module.
"""
warnings.warn('Input variables will be reordered in 1.0 to be (heights, u, v, f, dx, dy).'
'To update to new input format before 1.0 is released, use'
'`from metpy.future import ageostrophic_wind`.', FutureWarning)
from ..future import ageostrophic_wind as _ageostrophic_wind
return _ageostrophic_wind(heights, u, v, f, dx, dy, dim_order=dim_order)
u_geostrophic, v_geostrophic = geostrophic_wind(heights, f, dx, dy, dim_order=dim_order)
return u - u_geostrophic, v - v_geostrophic


@exporter.export
Expand Down Expand Up @@ -534,9 +533,9 @@ def montgomery_streamfunction(height, temperature):

@exporter.export
@preprocess_xarray
@check_units('[speed]', '[speed]', '[length]', '[length]', '[length]',
'[speed]', '[speed]')
def storm_relative_helicity(u, v, heights, depth, bottom=0 * units.m,
@check_units('[length]', '[speed]', '[speed]', '[length]',
bottom='[length]', storm_u='[speed]', storm_v='[speed]')
def storm_relative_helicity(heights, u, v, depth, *, bottom=0 * units.m,
storm_u=0 * units('m/s'), storm_v=0 * units('m/s')):
# Partially adapted from similar SharpPy code
r"""Calculate storm relative helicity.
Expand Down Expand Up @@ -577,13 +576,26 @@ def storm_relative_helicity(u, v, heights, depth, bottom=0 * units.m,
total storm-relative helicity
"""
warnings.warn('Input variables will be reordered in 1.0 to be (heights, u, v, depth, '
'bottom, storm_u, storm_v). To update to new input format before 1.0 is '
'released, use `from metpy.future import storm_relative_helicity`.',
FutureWarning)
from ..future import storm_relative_helicity as _storm_relative_helicity
return _storm_relative_helicity(heights, u, v, depth,
bottom=bottom, storm_u=storm_u, storm_v=storm_v)
_, u, v = get_layer_heights(heights, depth, u, v, with_agl=True, bottom=bottom)

storm_relative_u = u - storm_u
storm_relative_v = v - storm_v

int_layers = (storm_relative_u[1:] * storm_relative_v[:-1]
- storm_relative_u[:-1] * storm_relative_v[1:])

# Need to manually check for masked value because sum() on masked array with non-default
# mask will return a masked value rather than 0. See numpy/numpy#11736
positive_srh = int_layers[int_layers.magnitude > 0.].sum()
if np.ma.is_masked(positive_srh):
positive_srh = 0.0 * units('meter**2 / second**2')
negative_srh = int_layers[int_layers.magnitude < 0.].sum()
if np.ma.is_masked(negative_srh):
negative_srh = 0.0 * units('meter**2 / second**2')

return (positive_srh.to('meter ** 2 / second ** 2'),
negative_srh.to('meter ** 2 / second ** 2'),
(positive_srh + negative_srh).to('meter ** 2 / second ** 2'))


@exporter.export
Expand Down
182 changes: 1 addition & 181 deletions src/metpy/future.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,184 +3,4 @@
# SPDX-License-Identifier: BSD-3-Clause
"""Contains updated functions that will be modified in a future release."""

import numpy as np

from . import constants as mpconsts
from .calc.kinematics import ensure_yx_order, geostrophic_wind
from .calc.thermo import mixing_ratio, saturation_vapor_pressure
from .calc.tools import _remove_nans, get_layer, get_layer_heights
from .package_tools import Exporter
from .units import check_units, units
from .xarray import preprocess_xarray


exporter = Exporter(globals())


@exporter.export
@preprocess_xarray
@ensure_yx_order
@check_units(f='[frequency]', u='[speed]', v='[speed]', dx='[length]', dy='[length]')
def ageostrophic_wind(heights, u, v, f, dx, dy, dim_order='yx'):
r"""Calculate the ageostrophic wind given from the heights or geopotential.
Parameters
----------
heights : (M, N) ndarray
The height or geopotential field.
u : (M, N) `pint.Quantity`
The u wind field.
v : (M, N) `pint.Quantity`
The u wind field.
f : array_like
The coriolis parameter. This can be a scalar to be applied
everywhere or an array of values.
dx : `pint.Quantity`
The grid spacing(s) in the x-direction. If an array, there should be one item less than
the size of `heights` along the applicable axis.
dy : `pint.Quantity`
The grid spacing(s) in the y-direction. If an array, there should be one item less than
the size of `heights` along the applicable axis.
Returns
-------
A 2-item tuple of arrays
A tuple of the u-component and v-component of the ageostrophic wind.
Notes
-----
If inputs have more than two dimensions, they are assumed to have either leading dimensions
of (x, y) or trailing dimensions of (y, x), depending on the value of ``dim_order``.
This function contains an updated input variable order from the same function in the
kinematics module. This version will be fully implemented in 1.0 and moved from the
`future` module back to the `kinematics` module.
"""
u_geostrophic, v_geostrophic = geostrophic_wind(heights, f, dx, dy, dim_order=dim_order)
return u - u_geostrophic, v - v_geostrophic


@exporter.export
@preprocess_xarray
@check_units('[pressure]', '[temperature]', bottom='[pressure]', top='[pressure]')
def precipitable_water(pressure, dewpt, *, bottom=None, top=None):
r"""Calculate precipitable water through the depth of a sounding.
Formula used is:
.. math:: -\frac{1}{\rho_l g} \int\limits_{p_\text{bottom}}^{p_\text{top}} r dp
from [Salby1996]_, p. 28.
Parameters
----------
pressure : `pint.Quantity`
Atmospheric pressure profile
dewpt : `pint.Quantity`
Atmospheric dewpoint profile
bottom: `pint.Quantity`, optional
Bottom of the layer, specified in pressure. Defaults to None (highest pressure).
top: `pint.Quantity`, optional
The top of the layer, specified in pressure. Defaults to None (lowest pressure).
Returns
-------
`pint.Quantity`
The precipitable water in the layer
Examples
--------
>>> pressure = np.array([1000, 950, 900]) * units.hPa
>>> dewpoint = np.array([20, 15, 10]) * units.degC
>>> pw = precipitable_water(pressure, dewpoint)
"""
# Sort pressure and dewpoint to be in decreasing pressure order (increasing height)
sort_inds = np.argsort(pressure)[::-1]
pressure = pressure[sort_inds]
dewpt = dewpt[sort_inds]

pressure, dewpt = _remove_nans(pressure, dewpt)

if top is None:
top = np.nanmin(pressure.magnitude) * pressure.units

if bottom is None:
bottom = np.nanmax(pressure.magnitude) * pressure.units

pres_layer, dewpt_layer = get_layer(pressure, dewpt, bottom=bottom, depth=bottom - top)

w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)

# Since pressure is in decreasing order, pw will be the opposite sign of that expected.
pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units)
/ (mpconsts.g * mpconsts.rho_l))
return pw.to('millimeters')


@exporter.export
@preprocess_xarray
@check_units('[length]', '[speed]', '[speed]', '[length]',
bottom='[length]', storm_u='[speed]', storm_v='[speed]')
def storm_relative_helicity(heights, u, v, depth, *, bottom=0 * units.m,
storm_u=0 * units('m/s'), storm_v=0 * units('m/s')):
# Partially adapted from similar SharpPy code
r"""Calculate storm relative helicity.
Calculates storm relatively helicity following [Markowski2010]_ 230-231.
.. math:: \int\limits_0^d (\bar v - c) \cdot \bar\omega_{h} \,dz
This is applied to the data from a hodograph with the following summation:
.. math:: \sum_{n = 1}^{N-1} [(u_{n+1} - c_{x})(v_{n} - c_{y}) -
(u_{n} - c_{x})(v_{n+1} - c_{y})]
Parameters
----------
u : array-like
u component winds
v : array-like
v component winds
heights : array-like
atmospheric heights, will be converted to AGL
depth : number
depth of the layer
bottom : number
height of layer bottom AGL (default is surface)
storm_u : number
u component of storm motion (default is 0 m/s)
storm_v : number
v component of storm motion (default is 0 m/s)
Returns
-------
`pint.Quantity`
positive storm-relative helicity
`pint.Quantity`
negative storm-relative helicity
`pint.Quantity`
total storm-relative helicity
"""
_, u, v = get_layer_heights(heights, depth, u, v, with_agl=True, bottom=bottom)

storm_relative_u = u - storm_u
storm_relative_v = v - storm_v

int_layers = (storm_relative_u[1:] * storm_relative_v[:-1]
- storm_relative_u[:-1] * storm_relative_v[1:])

# Need to manually check for masked value because sum() on masked array with non-default
# mask will return a masked value rather than 0. See numpy/numpy#11736
positive_srh = int_layers[int_layers.magnitude > 0.].sum()
if np.ma.is_masked(positive_srh):
positive_srh = 0.0 * units('meter**2 / second**2')
negative_srh = int_layers[int_layers.magnitude < 0.].sum()
if np.ma.is_masked(negative_srh):
negative_srh = 0.0 * units('meter**2 / second**2')

return (positive_srh.to('meter ** 2 / second ** 2'),
negative_srh.to('meter ** 2 / second ** 2'),
(positive_srh + negative_srh).to('meter ** 2 / second ** 2'))
# This file left as a placeholder for function changes in the future.

0 comments on commit 67eeec7

Please sign in to comment.