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-04 20:30 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	57
rs	9.6818

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

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

    return u_shr, v_shr


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

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