supercell_composite() - Code Metrics - Inspection of "Supercell Composite and Significant Tornado Parame..." - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#545)

unknown

created 2017-08-23 18:05 UTC

supercell_composite() B

↳ 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	40
rs	8.8571

# 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 [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, esrh, ebwd):
    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) * (esrh / 50 m^2/s^2) * (ebwd / 20 m/s)

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

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

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

    """
    ebwd = ebwd.to('m/s')
    ebwd = ebwd.magnitude
    inds = np.where((ebwd <= 20.) & (ebwd >= 10.))
    inds1 = np.where(ebwd < 10.)
    inds2 = np.where(ebwd > 20.)
    ebwd[inds] = ebwd[inds] / 20.
    ebwd[inds1] = 0.
    ebwd[inds2] = 1.
    sup_comp = (mucape.magnitude / 1000.) * (esrh.magnitude / 50.) * ebwd

    return sup_comp


@exporter.export
def significant_tornado(sbcape, sblcl, srh1, shr6):
    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) * (srh1 / 150 m^s/s^2) *
             (shr6 / 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
    srh1 : `pint.Quantity`
        Surface-1km storm-relative helicity
    shr6 : `pint.Quantity`
        Surface-6km bulk shear

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

    """
    shr6 = shr6.to('m/s')
    shr6 = shr6.magnitude
    sblcl = sblcl.magnitude
    ind = np.where((sblcl <= 2000.) & (sblcl >= 1000.))
    ind1 = np.where(sblcl < 1000.)
    ind2 = np.where(sblcl > 2000.)
    sind = np.where((shr6 <= 30.) & (shr6 >= 12.5))
    sind1 = np.where(shr6 < 12.5)
    sind2 = np.where(shr6 > 30.)
    sblcl[ind] = (2000. - sblcl[ind]) / 1000.
    sblcl[ind1] = 1.
    sblcl[ind2] = 0.
    shr6[sind] = shr6[sind] / 20.
    shr6[sind1] = 0.
    shr6[sind2] = 1.5
    sigtor = (sbcape.magnitude / 1500.) * sblcl * (srh1.magnitude / 150.) * shr6

    return sigtor


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 [Bunkers2000]_.
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
194			v_shr: `pint.Quantity`
195			v-component of layer bulk shear
196
197			"""
198			_, u_layer, v_layer = get_layer(pressure, u, v, heights=heights,
199			bottom=bottom, depth=depth)
200
201			u_shr = u_layer[-1] - u_layer[0]
202			v_shr = v_layer[-1] - v_layer[0]
203
204			return u_shr, v_shr
205
206
207			@exporter.export
208			def supercell_composite(mucape, esrh, ebwd):
209			r"""Calculate the supercell composite parameter.
210
211			The supercell composite parameter is designed to identify
212			environments favorable for the development of supercells,
213			and is calculated using the formula developed by
214			[Thompson2004]:
215
216			SCP = (mucape / 1000 J/kg) * (esrh / 50 m^2/s^2) * (ebwd / 20 m/s)
217
218			The ebwd term is set to zero below 10 m/s and
219			capped at 1 when ebwd exceeds 20 m/s.
220
221			Parameters
222			----------
223			mucape : `pint.Quantity`
224			Most-unstable CAPE
225			esrh : `pint.Quantity`
226			Effective-layer storm-relative helicity
227			ebwd : `pint.Quantity`
228			Effective bulk shear
229
230			Returns
231			-------
232			array-like
233			supercell composite
234
235			"""
236			ebwd = ebwd.to('m/s')
237			ebwd = ebwd.magnitude
238			inds = np.where((ebwd <= 20.) & (ebwd >= 10.))
239			inds1 = np.where(ebwd < 10.)
240			inds2 = np.where(ebwd > 20.)
241			ebwd[inds] = ebwd[inds] / 20.
242			ebwd[inds1] = 0.
243			ebwd[inds2] = 1.
244			sup_comp = (mucape.magnitude / 1000.) * (esrh.magnitude / 50.) * ebwd
245
246			return sup_comp
247
248
249			@exporter.export
250			def significant_tornado(sbcape, sblcl, srh1, shr6):
251			r"""Calculate the significant tornado parameter (fixed layer).
252
253			The significant tornado parameter is designed to identify
254			environments favorable for the production of significant
255			tornadoes contingent upon the development of supercells.
256			It's calculated according to the formula used on the SPC
257			mesoanalysis page, updated in [Thompson2004]:
258
259			sigtor = (sbcape / 1500 J/kg) * ((2000 m - sblcl) / 1000 m) * (srh1 / 150 m^s/s^2) *
260			(shr6 / 20 m/s)
261
262			The sblcl term is set to zero when the lcl is above 2000m and
263			capped at 1 when below 1000m, and the shr6 term is set to 0
264			when shr6 is below 12.5 m/s and maxed out at 1.5 when shr6
265			exceeds 30 m/s.
266
267			Parameters
268			----------
269			sbcape : `pint.Quantity`
270			Surface-based CAPE
271			sblcl : `pint.Quantity`
272			Surface-based lifted condensation level
273			srh1 : `pint.Quantity`
274			Surface-1km storm-relative helicity
275			shr6 : `pint.Quantity`
276			Surface-6km bulk shear
277
278			Returns
279			-------
280			array-like
281			significant tornado parameter
282
283			"""
284			shr6 = shr6.to('m/s')
285			shr6 = shr6.magnitude
286			sblcl = sblcl.magnitude
287			ind = np.where((sblcl <= 2000.) & (sblcl >= 1000.))
288			ind1 = np.where(sblcl < 1000.)
289			ind2 = np.where(sblcl > 2000.)
290			sind = np.where((shr6 <= 30.) & (shr6 >= 12.5))
291			sind1 = np.where(shr6 < 12.5)
292			sind2 = np.where(shr6 > 30.)
293			sblcl[ind] = (2000. - sblcl[ind]) / 1000.
294			sblcl[ind1] = 1.
295			sblcl[ind2] = 0.
296			shr6[sind] = shr6[sind] / 20.
297			shr6[sind1] = 0.
298			shr6[sind2] = 1.5
299			sigtor = (sbcape.magnitude / 1500.) * sblcl * (srh1.magnitude / 150.) * shr6
300
301			return sigtor
302

Unidata / MetPy

Pull Request — master (#545)

supercell_composite() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like