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 21:43 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	59
rs	9.597

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, concatenate, 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
    -------
    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

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

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

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

    # mean wind from 5.5-6km
    wind_5500m = concatenate(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
    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') / 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`
        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
    -------
    u_shr: `pint.Quantity'
        u-component of layer bulk shear, in m/s
    v_shr: `pint.Quantity'
        v-component of layer bulk shear, in m/s
    shr_mag: `pint.Quantity'
        magnitude of layer bulk shear, in m/s

    """
    _, 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]

    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, concatenate, 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			right_mover: `pint.Quantity`
127			U and v component of Bunkers RM storm motion
128			left_mover: `pint.Quantity`
129			U and v component of Bunkers LM storm motion
130			wind_mean: `pint.Quantity`
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			wind_mean = concatenate(mean_pressure_weighted(pressure, ums, vms, heights=heights,
139			depth=6000 * units('meter')))
140
141			# mean wind from sfc-500m
142			wind_500m = concatenate(mean_pressure_weighted(pressure, ums, vms, heights=heights,
143			depth=500 * units('meter')))
144
145			# mean wind from 5.5-6km
146			wind_5500m = concatenate(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			shear = wind_5500m - wind_500m
152
153			# Take the cross product of the wind shear and k, and divide by the vector magnitude and
154			# multiply by the deviaton empirically calculated in Bunkers (2000) (7.5 m/s)
155			shear_cross = concatenate([shear[1], -shear[0]])
156			rdev = shear_cross * (7.5 * units('m/s') / np.hypot(*shear))
157
158			# Add the deviations to the layer average wind to get the RM motion
159			right_mover = wind_mean + rdev
160
161			# Subtract the deviations to get the LM motion
162			left_mover = wind_mean - rdev
163
164			return right_mover, left_mover, wind_mean
165
166
167			@exporter.export
168			@check_units('[pressure]', '[speed]', '[speed]')
169			def bulk_shear(pressure, u, v, heights=None, bottom=None, depth=None):
170			r"""Calculate bulk shear through a layer.
171
172			Layer top and bottom specified in meters or pressure.
173
174			Parameters
175			----------
176			pressure : `pint.Quantity`
177			Atmospheric pressure profile
178			u : `pint.Quantity`
179			U-component of wind.
180			v : `pint.Quantity`
181			V-component of wind.
182			height : `pint.Quantity`
183			Heights from sounding
184			depth: `pint.Quantity`
185			The depth of the layer in meters or hPa
186			bottom: `pint.Quantity`
187			The bottom of the layer in meters or hPa.
188			If in meters, must be in the same coordinates as the given
189			heights (i.e., don't use meters AGL unless given heights
190			are in meters AGL.) Default is the surface (1st observation.)
191
192			Returns
193			-------
194			u_shr: `pint.Quantity'
195			u-component of layer bulk shear, in m/s
196			v_shr: `pint.Quantity'
197			v-component of layer bulk shear, in m/s
198			shr_mag: `pint.Quantity'
199			magnitude of layer bulk shear, in m/s
200
201			"""
202			_, u_layer, v_layer = get_layer(pressure, u, v, heights=heights,
203			bottom=bottom, depth=depth)
204
205			u_shr = u_layer[-1] - u_layer[0]
206			v_shr = v_layer[-1] - v_layer[0]
207
208			shr_mag = np.sqrt((u_shr 2) + (v_shr 2))
209
210			return u_shr, v_shr, shr_mag
211

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