bunkers_storm_motion() - Code Metrics - Inspection of "Storm Motion and Shear" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#536)

unknown

created 2017-08-03 20:17 UTC

bunkers_storm_motion() A

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	1
c	1
b	0
f	0
dl	0
loc	73
rs	9.0675

How to fix Long Method

# Copyright (c) 2008-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, units

exporter = Exporter(globals())


@exporter.export
@check_units('[temperature]', '[pressure]', '[pressure]')
def precipitable_water(dewpt, pressure, top=400 * units('hPa')):
    r"""Calculate precipitable water through the depth of a sounding.

    Default layer depth is sfc-400 hPa. Formula used is:

    .. math:: \frac{1}{pg} \int\limits_0^d x \,dp

    from [Tsonis2008]_, p. 170.

    Parameters
    ----------
    dewpt : `pint.Quantity`
        Atmospheric dewpoint profile
    pressure : `pint.Quantity`
        Atmospheric pressure profile
    top: `pint.Quantity`, optional
        The top of the layer, specified in pressure. Defaults to 400 hPa.

    Returns
    -------
    `pint.Quantity`
        The precipitable water in the layer

    """
    sort_inds = np.argsort(pressure[::-1])
    pressure = pressure[sort_inds]
    dewpt = dewpt[sort_inds]

    pres_layer, dewpt_layer = get_layer(pressure, dewpt, depth=pressure[0] - top)

    w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)
    # Since pressure is in decreasing order, pw will be the negative of what we want.
    # Thus the *-1
    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 [Bunkers et al, 2000]_.

    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
    -------
    RM_vector
        U and v component of Bunkers RM storm motion
    LM_vector
        U and v component of Bunkers LM storm motion
    mean_vector
        U and v component of sfc-6km mean flow

    """
    ums = u.to('m/s')
    vms = v.to('m/s')

    # mean wind from sfc-6km
    u6m, v6m = mean_pressure_weighted(pressure, ums, vms, heights=heights,
                                      depth=6000 * units('meter'))

    # mean wind from sfc-500m
    u5m, v5m = mean_pressure_weighted(pressure, ums, vms, heights=heights,
                                      depth=500 * units('meter'))

    # mean wind from 5.5-6km
    u55m, v55m = mean_pressure_weighted(pressure, ums, vms, heights=heights,
                                        depth=500 * units('meter'),
                                        bottom=heights[0] + 5500 * units('meter'))

    # Calculate the shear vector from sfc-500m to 5.5-6km
    u6shr = u55m - u5m
    v6shr = v55m - v5m

    # making the shear vector
    vl = [u6shr.magnitude, v6shr.magnitude, 0]

    # Create a k hat vector
    vk = [0, 0, 1]

    # 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)
    rdev = np.cross(vl, vk) * (7.5 / (np.sqrt(u6shr.magnitude ** 2 + v6shr.magnitude ** 2)))

    # Add the deviations to the layer average wind to get the RM motion
    urm = u6m.magnitude + rdev[0]
    vrm = v6m.magnitude + rdev[1]

    # Subtract the deviations to get the LM motion
    ulm = u6m.magnitude - rdev[0]
    vlm = v6m.magnitude - rdev[1]

    rm_vector = np.asarray([urm, vrm]) * units('m/s').to(u.units)

    lm_vector = np.asarray([ulm, vlm]) * units('m/s').to(u.units)

    mean_vector = np.asarray([u6m.magnitude, v6m.magnitude]) * units('m/s').to(u.units)

    return rm_vector, lm_vector, mean_vector


@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`
        Heights from sounding
    depth: `pint.Quantity`
        The depth of the layer in meters or hPa
    bottom: `pint.Quantity`
        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
    -------
    `pint.Quantity'
        u_shr: u-component of layer bulk shear, in m/s
    `pint.Quantity'
        v_shr: v-component of layer bulk shear, in m/s
    `pint.Quantity'
        shr_mag: magnitude of layer bulk shear, in m/s

    """
    u = u.to('meters/second')
    v = v.to('meters/second')

    sort_inds = np.argsort(pressure[::-1])
    pressure = pressure[sort_inds]
    heights = heights[sort_inds]
    u = u[sort_inds]
    v = v[sort_inds]

    w_int = get_layer(pressure, u, v, heights=heights, bottom=bottom, depth=depth)

    u_shr = w_int[1][-1] - w_int[1][0]
    v_shr = w_int[2][-1] - w_int[2][0]

    shr_mag = np.sqrt((u_shr ** 2) + (v_shr ** 2))

    return u_shr, v_shr, shr_mag


1			# Copyright (c) 2008-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, units
12
13			exporter = Exporter(globals())
14
15
16			@exporter.export
17			@check_units('[temperature]', '[pressure]', '[pressure]')
18			def precipitable_water(dewpt, pressure, top=400 * units('hPa')):
19			r"""Calculate precipitable water through the depth of a sounding.
20
21			Default layer depth is sfc-400 hPa. Formula used is:
22
23			.. math:: \frac{1}{pg} \int\limits_0^d x \,dp
24
25			from [Tsonis2008]_, p. 170.
26
27			Parameters
28			----------
29			dewpt : `pint.Quantity`
30			Atmospheric dewpoint profile
31			pressure : `pint.Quantity`
32			Atmospheric pressure profile
33			top: `pint.Quantity`, optional
34			The top of the layer, specified in pressure. Defaults to 400 hPa.
35
36			Returns
37			-------
38			`pint.Quantity`
39			The precipitable water in the layer
40
41			"""
42			sort_inds = np.argsort(pressure[::-1])
43			pressure = pressure[sort_inds]
44			dewpt = dewpt[sort_inds]
45
46			pres_layer, dewpt_layer = get_layer(pressure, dewpt, depth=pressure[0] - top)
47
48			w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)
49			# Since pressure is in decreasing order, pw will be the negative of what we want.
50			# Thus the *-1
51			pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units) /
52			(g * rho_l))
53			return pw.to('millimeters')
54
55
56			@exporter.export
57			@check_units('[pressure]')
58			def mean_pressure_weighted(pressure, args, *kwargs):
59			r"""Calculate pressure-weighted mean of an arbitrary variable through a layer.
60
61			Layer top and bottom specified in height or pressure.
62
63			Parameters
64			----------
65			pressure : `pint.Quantity`
66			Atmospheric pressure profile
67			*args : `pint.Quantity`
68			Parameters for which the pressure-weighted mean is to be calculated.
69			heights : `pint.Quantity`, optional
70			Heights from sounding. Standard atmosphere heights assumed (if needed)
71			if no heights are given.
72			bottom: `pint.Quantity`, optional
73			The bottom of the layer in either the provided height coordinate
74			or in pressure. Don't provide in meters AGL unless the provided
75			height coordinate is meters AGL. Default is the first observation,
76			assumed to be the surface.
77			depth: `pint.Quantity`, optional
78			The depth of the layer in meters or hPa.
79
80			Returns
81			-------
82			`pint.Quantity`
83			u_mean: u-component of layer mean wind.
84			`pint.Quantity`
85			v_mean: v-component of layer mean wind.
86
87			"""
88			heights = kwargs.pop('heights', None)
89			bottom = kwargs.pop('bottom', None)
90			depth = kwargs.pop('depth', None)
91			ret = [] # Returned variable means in layer
92			layer_arg = get_layer(pressure, *args, heights=heights,
93			bottom=bottom, depth=depth)
94			layer_p = layer_arg[0]
95			layer_arg = layer_arg[1:]
96			# Taking the integral of the weights (pressure) to feed into the weighting
97			# function. Said integral works out to this function:
98			pres_int = 0.5 * (layer_p[-1].magnitude2 - layer_p[0].magnitude2)
99			for i, datavar in enumerate(args):
100			arg_mean = np.trapz(layer_arg[i] * layer_p, x=layer_p) / pres_int
101			ret.append(arg_mean * datavar.units)
102
103			return ret
104
105
106			@exporter.export
107			@check_units('[pressure]', '[speed]', '[speed]', '[length]')
108			def bunkers_storm_motion(pressure, u, v, heights):
109			r"""Calculate the Bunkers right-mover and left-mover storm motions and sfc-6km mean flow.
110
111			Uses the storm motion calculation from [Bunkers et al, 2000]_.
112
113			Parameters
114			----------
115			pressure : array-like
116			Pressure from sounding
117			u : array-like
118			U component of the wind
119			v : array-like
120			V component of the wind
121			heights : array-like
122			Heights from sounding
123
124			Returns
125			-------
126			RM_vector
127			U and v component of Bunkers RM storm motion
128			LM_vector
129			U and v component of Bunkers LM storm motion
130			mean_vector
131			U and v component of sfc-6km mean flow
132
133			"""
134			ums = u.to('m/s')
135			vms = v.to('m/s')
136
137			# mean wind from sfc-6km
138			u6m, v6m = mean_pressure_weighted(pressure, ums, vms, heights=heights,
139			depth=6000 * units('meter'))
140
141			# mean wind from sfc-500m
142			u5m, v5m = mean_pressure_weighted(pressure, ums, vms, heights=heights,
143			depth=500 * units('meter'))
144
145			# mean wind from 5.5-6km
146			u55m, v55m = mean_pressure_weighted(pressure, ums, vms, heights=heights,
147			depth=500 * units('meter'),
148			bottom=heights[0] + 5500 * units('meter'))
149
150			# Calculate the shear vector from sfc-500m to 5.5-6km
151			u6shr = u55m - u5m
152			v6shr = v55m - v5m
153
154			# making the shear vector
155			vl = [u6shr.magnitude, v6shr.magnitude, 0]
156
157			# Create a k hat vector
158			vk = [0, 0, 1]
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			rdev = np.cross(vl, vk) * (7.5 / (np.sqrt(u6shr.magnitude 2 + v6shr.magnitude 2)))
163
164			# Add the deviations to the layer average wind to get the RM motion
165			urm = u6m.magnitude + rdev[0]
166			vrm = v6m.magnitude + rdev[1]
167
168			# Subtract the deviations to get the LM motion
169			ulm = u6m.magnitude - rdev[0]
170			vlm = v6m.magnitude - rdev[1]
171
172			rm_vector = np.asarray([urm, vrm]) * units('m/s').to(u.units)
173
174			lm_vector = np.asarray([ulm, vlm]) * units('m/s').to(u.units)
175
176			mean_vector = np.asarray([u6m.magnitude, v6m.magnitude]) * units('m/s').to(u.units)
177
178			return rm_vector, lm_vector, mean_vector
179
180
181			@exporter.export
182			@check_units('[pressure]', '[speed]', '[speed]')
183			def bulk_shear(pressure, u, v, heights=None, bottom=None, depth=None):
184			r"""Calculate bulk shear through a layer.
185
186			Layer top and bottom specified in meters or pressure.
187
188			Parameters
189			----------
190			pressure : `pint.Quantity`
191			Atmospheric pressure profile
192			u : `pint.Quantity`
193			U-component of wind.
194			v : `pint.Quantity`
195			V-component of wind.
196			height : `pint.Quantity`
197			Heights from sounding
198			depth: `pint.Quantity`
199			The depth of the layer in meters or hPa
200			bottom: `pint.Quantity`
201			The bottom of the layer in meters or hPa.
202			If in meters, must be in the same coordinates as the given
203			heights (i.e., don't use meters AGL unless given heights
204			are in meters AGL.) Default is the surface (1st observation.)
205
206			Returns
207			-------
208			`pint.Quantity'
209			u_shr: u-component of layer bulk shear, in m/s
210			`pint.Quantity'
211			v_shr: v-component of layer bulk shear, in m/s
212			`pint.Quantity'
213			shr_mag: magnitude of layer bulk shear, in m/s
214
215			"""
216			u = u.to('meters/second')
217			v = v.to('meters/second')
218
219			sort_inds = np.argsort(pressure[::-1])
220			pressure = pressure[sort_inds]
221			heights = heights[sort_inds]
222			u = u[sort_inds]
223			v = v[sort_inds]
224
225			w_int = get_layer(pressure, u, v, heights=heights, bottom=bottom, depth=depth)
226
227			u_shr = w_int[1][-1] - w_int[1][0]
228			v_shr = w_int[2][-1] - w_int[2][0]
229
230			shr_mag = np.sqrt((u_shr 2) + (v_shr 2))
231
232			return u_shr, v_shr, shr_mag
233

Unidata / MetPy

Pull Request — master (#536)

bunkers_storm_motion() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like