precipitable_water() - Code Metrics - Inspection of "Merge pull request #595 from jrleeman/precip_water..." - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — master ( f1beb3...24e483 )

by Ryan

created 2017-10-30 20:31 UTC

precipitable_water() A

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	7
Ratio	14.89 %

Importance

Changes	1
Bugs	0	Features	1

Metric	Value
cc	3
c	1
b	0
f	1
dl	7
loc	47
rs	9.0303

# Copyright (c) 2017 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains calculation of various derived indices."""
import numpy as np

from .thermo import mixing_ratio, saturation_vapor_pressure
from .tools import get_layer
from ..constants import g, rho_l
from ..package_tools import Exporter
from ..units import check_units, concatenate, units

exporter = Exporter(globals())


@exporter.export
@check_units('[temperature]', '[pressure]', '[pressure]')
def precipitable_water(dewpt, pressure, 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
    ----------
    dewpt : `pint.Quantity`
        Atmospheric dewpoint profile
    pressure : `pint.Quantity`
        Atmospheric pressure 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

    """
    # 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]

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

    if bottom is None:
        bottom = np.nanmax(pressure) * 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) /
                (g * rho_l))
    return pw.to('millimeters')


@exporter.export
@check_units('[pressure]')
def mean_pressure_weighted(pressure, *args, **kwargs):
    r"""Calculate pressure-weighted mean of an arbitrary variable through a layer.

    Layer top and bottom specified in height or pressure.

    Parameters
    ----------
    pressure : `pint.Quantity`
        Atmospheric pressure profile
    *args : `pint.Quantity`
        Parameters for which the pressure-weighted mean is to be calculated.
    heights : `pint.Quantity`, optional
        Heights from sounding. Standard atmosphere heights assumed (if needed)
        if no heights are given.
    bottom: `pint.Quantity`, optional
        The bottom of the layer in either the provided height coordinate
        or in pressure. Don't provide in meters AGL unless the provided
        height coordinate is meters AGL. Default is the first observation,
        assumed to be the surface.
    depth: `pint.Quantity`, optional
        The depth of the layer in meters or hPa.

    Returns
    -------
    `pint.Quantity`
        u_mean: u-component of layer mean wind.
    `pint.Quantity`
        v_mean: v-component of layer mean wind.

    """
    heights = kwargs.pop('heights', None)
    bottom = kwargs.pop('bottom', None)
    depth = kwargs.pop('depth', None)
    ret = []  # Returned variable means in layer
    layer_arg = get_layer(pressure, *args, heights=heights,
                          bottom=bottom, depth=depth)
    layer_p = layer_arg[0]
    layer_arg = layer_arg[1:]
    # Taking the integral of the weights (pressure) to feed into the weighting
    # function. Said integral works out to this function:
    pres_int = 0.5 * (layer_p[-1].magnitude**2 - layer_p[0].magnitude**2)
    for i, datavar in enumerate(args):
        arg_mean = np.trapz(layer_arg[i] * layer_p, x=layer_p) / pres_int
        ret.append(arg_mean * datavar.units)

    return ret


@exporter.export
@check_units('[pressure]', '[speed]', '[speed]', '[length]')
def bunkers_storm_motion(pressure, u, v, heights):
    r"""Calculate the Bunkers right-mover and left-mover storm motions and sfc-6km mean flow.

    Uses the storm motion calculation from [Bunkers2000]_.

    Parameters
    ----------
    pressure : array-like
        Pressure from sounding
    u : array-like
        U component of the wind
    v : array-like
        V component of the wind
    heights : array-like
        Heights from sounding

    Returns
    -------
    right_mover: `pint.Quantity`
        U and v component of Bunkers RM storm motion
    left_mover: `pint.Quantity`
        U and v component of Bunkers LM storm motion
    wind_mean: `pint.Quantity`
        U and v component of sfc-6km mean flow

    """
    # mean wind from sfc-6km
    wind_mean = concatenate(mean_pressure_weighted(pressure, u, v, heights=heights,
                                                   depth=6000 * units('meter')))

    # mean wind from sfc-500m
    wind_500m = concatenate(mean_pressure_weighted(pressure, u, v, heights=heights,
                                                   depth=500 * units('meter')))

    # mean wind from 5.5-6km
    wind_5500m = concatenate(mean_pressure_weighted(pressure, u, v, heights=heights,
                                                    depth=500 * units('meter'),
                                                    bottom=heights[0] +
                                                    5500 * units('meter')))

    # Calculate the shear vector from sfc-500m to 5.5-6km
    shear = wind_5500m - wind_500m

    # Take the cross product of the wind shear and k, and divide by the vector magnitude and
    # multiply by the deviaton empirically calculated in Bunkers (2000) (7.5 m/s)
    shear_cross = concatenate([shear[1], -shear[0]])
    rdev = shear_cross * (7.5 * units('m/s').to(u.units) / np.hypot(*shear))

    # Add the deviations to the layer average wind to get the RM motion
    right_mover = wind_mean + rdev

    # Subtract the deviations to get the LM motion
    left_mover = wind_mean - rdev

    return right_mover, left_mover, wind_mean


@exporter.export
@check_units('[pressure]', '[speed]', '[speed]')
def bulk_shear(pressure, u, v, heights=None, bottom=None, depth=None):
    r"""Calculate bulk shear through a layer.

    Layer top and bottom specified in meters or pressure.

    Parameters
    ----------
    pressure : `pint.Quantity`
        Atmospheric pressure profile
    u : `pint.Quantity`
        U-component of wind.
    v : `pint.Quantity`
        V-component of wind.
    height : `pint.Quantity`, optional
        Heights from sounding
    depth: `pint.Quantity`, optional
        The depth of the layer in meters or hPa
    bottom: `pint.Quantity`, optional
        The bottom of the layer in meters or hPa.
        If in meters, must be in the same coordinates as the given
        heights (i.e., don't use meters AGL unless given heights
        are in meters AGL.) Default is the surface (1st observation.)

    Returns
    -------
    u_shr: `pint.Quantity`
        u-component of layer bulk shear
    v_shr: `pint.Quantity`
        v-component of layer bulk shear

    """
    _, u_layer, v_layer = get_layer(pressure, u, v, heights=heights,
                                    bottom=bottom, depth=depth)

    u_shr = u_layer[-1] - u_layer[0]
    v_shr = v_layer[-1] - v_layer[0]

    return u_shr, v_shr


@exporter.export
def supercell_composite(mucape, effective_storm_helicity, effective_shear):
    r"""Calculate the supercell composite parameter.

    The supercell composite parameter is designed to identify
    environments favorable for the development of supercells,
    and is calculated using the formula developed by
    [Thompson2004]_:

    SCP = (mucape / 1000 J/kg) * (effective_storm_helicity / 50 m^2/s^2) *
          (effective_shear / 20 m/s)

    The effective_shear term is set to zero below 10 m/s and
    capped at 1 when effective_shear exceeds 20 m/s.

    Parameters
    ----------
    mucape : `pint.Quantity`
        Most-unstable CAPE
    effective_storm_helicity : `pint.Quantity`
        Effective-layer storm-relative helicity
    effective_shear : `pint.Quantity`
        Effective bulk shear

    Returns
    -------
    array-like
        supercell composite

    """
    effective_shear = np.clip(effective_shear, None, 20 * units('m/s'))
    effective_shear[effective_shear < 10 * units('m/s')] = 0 * units('m/s')
    effective_shear = effective_shear / (20 * units('m/s'))

    return ((mucape / (1000 * units('J/kg'))) *
            (effective_storm_helicity / (50 * units('m^2/s^2'))) *
            effective_shear).to('dimensionless')


@exporter.export
def significant_tornado(sbcape, sblcl, storm_helicity_1km, shear_6km):
    r"""Calculate the significant tornado parameter (fixed layer).

    The significant tornado parameter is designed to identify
    environments favorable for the production of significant
    tornadoes contingent upon the development of supercells.
    It's calculated according to the formula used on the SPC
    mesoanalysis page, updated in [Thompson2004]_:

    sigtor = (sbcape / 1500 J/kg) * ((2000 m - sblcl) / 1000 m) *
             (storm_helicity_1km / 150 m^s/s^2) * (shear_6km6 / 20 m/s)

    The sblcl term is set to zero when the lcl is above 2000m and
    capped at 1 when below 1000m, and the shr6 term is set to 0
    when shr6 is below 12.5 m/s and maxed out at 1.5 when shr6
    exceeds 30 m/s.

    Parameters
    ----------
    sbcape : `pint.Quantity`
        Surface-based CAPE
    sblcl : `pint.Quantity`
        Surface-based lifted condensation level
    storm_helicity_1km : `pint.Quantity`
        Surface-1km storm-relative helicity
    shear_6km : `pint.Quantity`
        Surface-6km bulk shear

    Returns
    -------
    array-like
        significant tornado parameter

    """
    sblcl = np.clip(sblcl, 1000 * units('meter'), 2000 * units('meter'))
    sblcl[sblcl > 2000 * units('meter')] = 0 * units('meter')
    sblcl = (2000. * units('meter') - sblcl) / (1000. * units('meter'))
    shear_6km = np.clip(shear_6km, None, 30 * units('m/s'))
    shear_6km[shear_6km < 12.5 * units('m/s')] = 0 * units('m/s')
    shear_6km = shear_6km / (20 * units('m/s'))

    return ((sbcape / (1500. * units('J/kg'))) *
            sblcl * (storm_helicity_1km / (150. * units('m^2/s^2'))) * shear_6km)


1		# Copyright (c) 2017 MetPy Developers.
2		# Distributed under the terms of the BSD 3-Clause License.
3		# SPDX-License-Identifier: BSD-3-Clause
4		"""Contains calculation of various derived indices."""
5		import numpy as np
6
7		from .thermo import mixing_ratio, saturation_vapor_pressure
8		from .tools import get_layer
9		from ..constants import g, rho_l
10		from ..package_tools import Exporter
11		from ..units import check_units, concatenate, units
12
13		exporter = Exporter(globals())
14
15
16		@exporter.export
17		@check_units('[temperature]', '[pressure]', '[pressure]')
18		def precipitable_water(dewpt, pressure, bottom=None, top=None):
19		r"""Calculate precipitable water through the depth of a sounding.
20
21		Formula used is:
22
23		.. math:: -\frac{1}{\rho_l g} \int\limits_{p_\text{bottom}}^{p_\text{top}} r dp
24
25		from [Salby1996]_, p. 28.
26
27		Parameters
28		----------
29		dewpt : `pint.Quantity`
30		Atmospheric dewpoint profile
31		pressure : `pint.Quantity`
32		Atmospheric pressure profile
33		bottom: `pint.Quantity`, optional
34		Bottom of the layer, specified in pressure. Defaults to None (highest pressure).
35		top: `pint.Quantity`, optional
36		The top of the layer, specified in pressure. Defaults to None (lowest pressure).
37
38		Returns
39		-------
40		`pint.Quantity`
41		The precipitable water in the layer
42
43		"""
44		# Sort pressure and dewpoint to be in decreasing pressure order (increasing height)
45		sort_inds = np.argsort(pressure)[::-1]
46		pressure = pressure[sort_inds]
47		dewpt = dewpt[sort_inds]
48
49		if top is None:
50		top = np.nanmin(pressure) * pressure.units
51
52		if bottom is None:
53		bottom = np.nanmax(pressure) * pressure.units
54
55		pres_layer, dewpt_layer = get_layer(pressure, dewpt, bottom=bottom, depth=bottom - top)
56	View Code Duplication
		0 ignored issues – show Duplication introduced 2017-10-26 18:34 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
57		w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)
58
59		# Since pressure is in decreasing order, pw will be the opposite sign of that expected.
60		pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units) /
61		(g * rho_l))
62		return pw.to('millimeters')
63
64
65		@exporter.export
66		@check_units('[pressure]')
67		def mean_pressure_weighted(pressure, args, *kwargs):
68		r"""Calculate pressure-weighted mean of an arbitrary variable through a layer.
69
70		Layer top and bottom specified in height or pressure.
71
72		Parameters
73		----------
74		pressure : `pint.Quantity`
75		Atmospheric pressure profile
76		*args : `pint.Quantity`
77		Parameters for which the pressure-weighted mean is to be calculated.
78		heights : `pint.Quantity`, optional
79		Heights from sounding. Standard atmosphere heights assumed (if needed)
80		if no heights are given.
81		bottom: `pint.Quantity`, optional
82		The bottom of the layer in either the provided height coordinate
83		or in pressure. Don't provide in meters AGL unless the provided
84		height coordinate is meters AGL. Default is the first observation,
85		assumed to be the surface.
86		depth: `pint.Quantity`, optional
87		The depth of the layer in meters or hPa.
88
89		Returns
90		-------
91		`pint.Quantity`
92		u_mean: u-component of layer mean wind.
93		`pint.Quantity`
94		v_mean: v-component of layer mean wind.
95
96		"""
97		heights = kwargs.pop('heights', None)
98		bottom = kwargs.pop('bottom', None)
99		depth = kwargs.pop('depth', None)
100		ret = [] # Returned variable means in layer
101		layer_arg = get_layer(pressure, *args, heights=heights,
102		bottom=bottom, depth=depth)
103		layer_p = layer_arg[0]
104		layer_arg = layer_arg[1:]
105		# Taking the integral of the weights (pressure) to feed into the weighting
106		# function. Said integral works out to this function:
107		pres_int = 0.5 * (layer_p[-1].magnitude2 - layer_p[0].magnitude2)
108		for i, datavar in enumerate(args):
109		arg_mean = np.trapz(layer_arg[i] * layer_p, x=layer_p) / pres_int
110		ret.append(arg_mean * datavar.units)
111
112		return ret
113
114
115		@exporter.export
116		@check_units('[pressure]', '[speed]', '[speed]', '[length]')
117		def bunkers_storm_motion(pressure, u, v, heights):
118		r"""Calculate the Bunkers right-mover and left-mover storm motions and sfc-6km mean flow.
119
120		Uses the storm motion calculation from [Bunkers2000]_.
121
122		Parameters
123		----------
124		pressure : array-like
125		Pressure from sounding
126		u : array-like
127		U component of the wind
128		v : array-like
129		V component of the wind
130		heights : array-like
131		Heights from sounding
132
133		Returns
134		-------
135		right_mover: `pint.Quantity`
136		U and v component of Bunkers RM storm motion
137		left_mover: `pint.Quantity`
138		U and v component of Bunkers LM storm motion
139		wind_mean: `pint.Quantity`
140		U and v component of sfc-6km mean flow
141
142		"""
143		# mean wind from sfc-6km
144		wind_mean = concatenate(mean_pressure_weighted(pressure, u, v, heights=heights,
145		depth=6000 * units('meter')))
146
147		# mean wind from sfc-500m
148		wind_500m = concatenate(mean_pressure_weighted(pressure, u, v, heights=heights,
149		depth=500 * units('meter')))
150
151		# mean wind from 5.5-6km
152		wind_5500m = concatenate(mean_pressure_weighted(pressure, u, v, heights=heights,
153		depth=500 * units('meter'),
154		bottom=heights[0] +
155		5500 * units('meter')))
156
157		# Calculate the shear vector from sfc-500m to 5.5-6km
158		shear = wind_5500m - wind_500m
159
160		# Take the cross product of the wind shear and k, and divide by the vector magnitude and
161		# multiply by the deviaton empirically calculated in Bunkers (2000) (7.5 m/s)
162		shear_cross = concatenate([shear[1], -shear[0]])
163		rdev = shear_cross * (7.5 * units('m/s').to(u.units) / np.hypot(*shear))
164
165		# Add the deviations to the layer average wind to get the RM motion
166		right_mover = wind_mean + rdev
167
168		# Subtract the deviations to get the LM motion
169		left_mover = wind_mean - rdev
170
171		return right_mover, left_mover, wind_mean
172
173
174		@exporter.export
175		@check_units('[pressure]', '[speed]', '[speed]')
176		def bulk_shear(pressure, u, v, heights=None, bottom=None, depth=None):
177		r"""Calculate bulk shear through a layer.
178
179		Layer top and bottom specified in meters or pressure.
180
181		Parameters
182		----------
183		pressure : `pint.Quantity`
184		Atmospheric pressure profile
185		u : `pint.Quantity`
186		U-component of wind.
187		v : `pint.Quantity`
188		V-component of wind.
189		height : `pint.Quantity`, optional
190		Heights from sounding
191		depth: `pint.Quantity`, optional
192		The depth of the layer in meters or hPa
193		bottom: `pint.Quantity`, optional
194		The bottom of the layer in meters or hPa.
195		If in meters, must be in the same coordinates as the given
196		heights (i.e., don't use meters AGL unless given heights
197		are in meters AGL.) Default is the surface (1st observation.)
198
199		Returns
200		-------
201		u_shr: `pint.Quantity`
202		u-component of layer bulk shear
203		v_shr: `pint.Quantity`
204		v-component of layer bulk shear
205
206		"""
207		_, u_layer, v_layer = get_layer(pressure, u, v, heights=heights,
208		bottom=bottom, depth=depth)
209
210		u_shr = u_layer[-1] - u_layer[0]
211		v_shr = v_layer[-1] - v_layer[0]
212
213		return u_shr, v_shr
214
215
216		@exporter.export
217		def supercell_composite(mucape, effective_storm_helicity, effective_shear):
218		r"""Calculate the supercell composite parameter.
219
220		The supercell composite parameter is designed to identify
221		environments favorable for the development of supercells,
222		and is calculated using the formula developed by
223		[Thompson2004]_:
224
225		SCP = (mucape / 1000 J/kg) * (effective_storm_helicity / 50 m^2/s^2) *
226		(effective_shear / 20 m/s)
227
228		The effective_shear term is set to zero below 10 m/s and
229		capped at 1 when effective_shear exceeds 20 m/s.
230
231		Parameters
232		----------
233		mucape : `pint.Quantity`
234		Most-unstable CAPE
235		effective_storm_helicity : `pint.Quantity`
236		Effective-layer storm-relative helicity
237		effective_shear : `pint.Quantity`
238		Effective bulk shear
239
240		Returns
241		-------
242		array-like
243		supercell composite
244
245		"""
246		effective_shear = np.clip(effective_shear, None, 20 * units('m/s'))
247		effective_shear[effective_shear < 10 * units('m/s')] = 0 * units('m/s')
248		effective_shear = effective_shear / (20 * units('m/s'))
249
250		return ((mucape / (1000 * units('J/kg'))) *
251		(effective_storm_helicity / (50 * units('m^2/s^2'))) *
252		effective_shear).to('dimensionless')
253
254
255		@exporter.export
256		def significant_tornado(sbcape, sblcl, storm_helicity_1km, shear_6km):
257		r"""Calculate the significant tornado parameter (fixed layer).
258
259		The significant tornado parameter is designed to identify
260		environments favorable for the production of significant
261		tornadoes contingent upon the development of supercells.
262		It's calculated according to the formula used on the SPC
263		mesoanalysis page, updated in [Thompson2004]_:
264
265		sigtor = (sbcape / 1500 J/kg) * ((2000 m - sblcl) / 1000 m) *
266		(storm_helicity_1km / 150 m^s/s^2) * (shear_6km6 / 20 m/s)
267
268		The sblcl term is set to zero when the lcl is above 2000m and
269		capped at 1 when below 1000m, and the shr6 term is set to 0
270		when shr6 is below 12.5 m/s and maxed out at 1.5 when shr6
271		exceeds 30 m/s.
272
273		Parameters
274		----------
275		sbcape : `pint.Quantity`
276		Surface-based CAPE
277		sblcl : `pint.Quantity`
278		Surface-based lifted condensation level
279		storm_helicity_1km : `pint.Quantity`
280		Surface-1km storm-relative helicity
281		shear_6km : `pint.Quantity`
282		Surface-6km bulk shear
283
284		Returns
285		-------
286		array-like
287		significant tornado parameter
288
289		"""
290		sblcl = np.clip(sblcl, 1000 * units('meter'), 2000 * units('meter'))
291		sblcl[sblcl > 2000 * units('meter')] = 0 * units('meter')
292		sblcl = (2000. * units('meter') - sblcl) / (1000. * units('meter'))
293		shear_6km = np.clip(shear_6km, None, 30 * units('m/s'))
294		shear_6km[shear_6km < 12.5 * units('m/s')] = 0 * units('m/s')
295		shear_6km = shear_6km / (20 * units('m/s'))
296
297		return ((sbcape / (1500. * units('J/kg'))) *
298		sblcl * (storm_helicity_1km / (150. * units('m^2/s^2'))) * shear_6km)
299

Unidata / MetPy

Push — master ( f1beb3...24e483 )

precipitable_water() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like