From 6c0c65e85c39e46e5b874f5b5dce1e31a02818d1 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Tue, 3 Oct 2017 07:51:38 -0600 Subject: [PATCH 1/7] [MNT] Cleanup storm_relative_helicity --- metpy/calc/kinematics.py | 56 ++++++++++++++++------------------------ 1 file changed, 22 insertions(+), 34 deletions(-) diff --git a/metpy/calc/kinematics.py b/metpy/calc/kinematics.py index 8310a185d76..0f68f05b5e4 100644 --- a/metpy/calc/kinematics.py +++ b/metpy/calc/kinematics.py @@ -541,16 +541,12 @@ def montgomery_streamfunction(height, temperature): @exporter.export @check_units('[speed]', '[speed]', '[pressure]', '[length]', '[length]', '[length]', '[speed]', '[speed]') -def storm_relative_helicity(u, v, p, hgt, top, bottom=0 * units('meter'), +def storm_relative_helicity(u, v, p, heights, top, bottom=0 * units('meter'), storm_u=0 * units('m/s'), storm_v=0 * units('m/s')): # Partially adapted from similar SharpPy code r"""Calculate storm relative helicity. - Needs u and v wind components, heights and pressures, - and top and bottom of SRH layer. An optional storm - motion vector can be specified. SRH is calculated using the - equation specified on p. 230-231 in the Markowski and Richardson - meso textbook [Markowski2010]. + Calculates storm relatively helicity following [Markowski2010] 230-231. .. math:: \int\limits_0^d (\bar v - c) \cdot \bar\omega_{h} \,dz @@ -562,18 +558,17 @@ def storm_relative_helicity(u, v, p, hgt, top, bottom=0 * units('meter'), Parameters ---------- u : array-like - The u components of winds, same length as hgts + u component winds v : array-like - The u components of winds, same length as hgts + v component winds p : array-like - Pressure in hPa, same length as hgts - hgt : array-like - The heights associatd with the data, provided in meters above mean - sea level and converted into meters AGL. + atmospheric pressure + heights : array-like + atmospheric heights top : number - The height of the top of the desired layer for SRH. + layer top bottom : number - The height at the bottom of the SRH layer. Default is sfc (None). + layer bottom (default is surface) storm_u : number u component of storm motion storm_v : number @@ -581,29 +576,22 @@ def storm_relative_helicity(u, v, p, hgt, top, bottom=0 * units('meter'), Returns ------- - number - p_srh : positive storm-relative helicity - number - n_srh : negative storm-relative helicity - number - t_srh : total storm-relative helicity + `pint.Quantity, pint.Quantity, pint.Quantity` + positive, negative, total storm-relative helicity """ - # Converting to m/s to make sure output is in m^2/s^2 - u = u.to('meters/second') - v = v.to('meters/second') - storm_u = storm_u.to('meters/second') - storm_v = storm_v.to('meters/second') + _, layer_u, layer_v = get_layer(p, u, v, heights=heights, + bottom=bottom, depth=top - bottom) - w_int = get_layer(p, u, v, heights=hgt, bottom=bottom, depth=top - bottom) + storm_relative_u = layer_u - storm_u + storm_relative_v = layer_v - storm_v - sru = w_int[1] - storm_u - srv = w_int[2] - storm_v + int_layers = (storm_relative_u[1:] * storm_relative_v[:-1] - + storm_relative_u[:-1] * storm_relative_v[1:]) - int_layers = sru[1:] * srv[:-1] - sru[:-1] * srv[1:] + positive_srh = int_layers[int_layers.magnitude > 0.].sum() + negative_srh = int_layers[int_layers.magnitude < 0.].sum() - p_srh = int_layers[int_layers.magnitude > 0.].sum() - n_srh = int_layers[int_layers.magnitude < 0.].sum() - t_srh = p_srh + n_srh - - return p_srh, n_srh, t_srh + return (positive_srh.to('meter ** 2 / second ** 2'), + negative_srh.to('meter ** 2 / second ** 2'), + (positive_srh + negative_srh).to('meter ** 2 / second ** 2')) From 612a2dc29b5a07b298becc70e9b957cf45632ca6 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Tue, 3 Oct 2017 13:08:21 -0600 Subject: [PATCH 2/7] Add get_layer_heights function for retreiving layers without pressure data --- metpy/calc/tools.py | 91 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/metpy/calc/tools.py b/metpy/calc/tools.py index 5971e8f33fc..99adc72ef18 100644 --- a/metpy/calc/tools.py +++ b/metpy/calc/tools.py @@ -394,6 +394,97 @@ def _get_bound_pressure_height(pressure, bound, heights=None, interpolate=True): return bound_pressure, bound_height +@exporter.export +@check_units('[length]') +def get_layer_heights(heights, depth, *args, **kwargs): + """Return an atmospheric layer from upper air data with the requested bottom and depth. + + This function will subset an upper air dataset to contain only the specified layer using + the heights only. + + Parameters + ---------- + heights : array-like + Atmospheric heights + depth : `pint.Quantity` + The thickness of the layer + *args : array-like + Atmospheric variable(s) measured at the given pressures + bottom : `pint.Quantity`, optional + The bottom of the layer + interpolate : bool, optional + Interpolate the top and bottom points if they are not in the given data. Defaults + to True. + with_agl : bool, optional + Returns the heights as above ground level by subtracting the minimum height in the + provided heights. Defaults to False. + + Returns + ------- + `pint.Quantity, pint.Quantity` + The height and data variables of the layer + + """ + bottom = kwargs.pop('bottom', None) + interpolate = kwargs.pop('interpolate', True) + with_agl = kwargs.pop('with_agl', False) + + # Make sure pressure and datavars are the same length + for datavar in args: + if len(heights) != len(datavar): + raise ValueError('Height and data variables must have the same length.') + + # If we want things in AGL, subtract the minimum height from all height values + if with_agl: + sfc_height = np.min(heights) + heights -= sfc_height + + # If the bottom is not specified, make it the surface + if bottom is None: + bottom = heights[0] + + # Make heights and arguments base units + heights = heights.to_base_units() + bottom = bottom.to_base_units() + + # Calculate the top of the layer + top = bottom + depth + + ret = [] # returned data variables in layer + + # Ensure heights are sorted in ascending order + sort_inds = np.argsort(heights) + heights = heights[sort_inds] + + # Mask based on top and bottom + inds = (heights >= bottom) & (heights <= top) + heights_interp = heights[inds] + + # Interpolate heights at bounds if necessary and sort + if interpolate: + # If we don't have the bottom or top requested, append them + if top not in heights_interp: + heights_interp = np.sort(np.append(heights_interp, top)) * heights.units + if bottom not in heights_interp: + heights_interp = np.sort(np.append(heights_interp, bottom)) * heights.units + + ret.append(heights_interp) + + for datavar in args: + # Ensure that things are sorted in ascending order + datavar = datavar[sort_inds] + + if interpolate: + # Interpolate for the possibly missing bottom/top values + datavar_interp = interp(heights_interp, heights, datavar) + datavar = datavar_interp + else: + datavar = datavar[inds] + + ret.append(datavar) + return ret + + @exporter.export @check_units('[pressure]') def get_layer(pressure, *args, **kwargs): From c48306ff23f0cdf487d2625dcefbfb6a954b171f Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Tue, 3 Oct 2017 14:45:13 -0600 Subject: [PATCH 3/7] Add tests for get_layer_heights --- metpy/calc/tests/test_tools.py | 49 +++++++++++++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/metpy/calc/tests/test_tools.py b/metpy/calc/tests/test_tools.py index 342e31b44d9..6ab65efcaf2 100644 --- a/metpy/calc/tests/test_tools.py +++ b/metpy/calc/tests/test_tools.py @@ -7,7 +7,8 @@ import numpy.ma as ma import pytest -from metpy.calc import (find_intersections, get_layer, interp, interpolate_nans, log_interp, +from metpy.calc import (find_intersections, get_layer, get_layer_heights, + interp, interpolate_nans, log_interp, nearest_intersection_idx, pressure_to_height_std, reduce_point_density, resample_nn_1d) from metpy.calc.tools import (_get_bound_pressure_height, _greater_or_close, _less_or_close, @@ -475,3 +476,49 @@ def test_less_or_close(): truth = np.array([True, True, True, True, True, False]) res = _less_or_close(x, comparison_value) assert_array_equal(res, truth) + + +def test_get_layer_heights_interpolation(): + """Test get_layer_heights with interpolation.""" + heights = np.arange(10) * units.km + data = heights.m * 2 * units.degC + heights, data = get_layer_heights(heights, 5000 * units.m, data, bottom=1500 * units.m) + heights_true = np.array([1.5, 2, 3, 4, 5, 6, 6.5]) * units.km + data_true = heights_true.m * 2 * units.degC + assert_array_almost_equal(heights_true, heights, 6) + assert_array_almost_equal(data_true, data, 6) + + +def test_get_layer_heights_no_interpolation(): + """Test get_layer_heights without interpolation.""" + heights = np.arange(10) * units.km + data = heights.m * 2 * units.degC + heights, data = get_layer_heights(heights, 5000 * units.m, data, + bottom=1500 * units.m, interpolate=False) + heights_true = np.array([2, 3, 4, 5, 6]) * units.km + data_true = heights_true.m * 2 * units.degC + assert_array_almost_equal(heights_true, heights, 6) + assert_array_almost_equal(data_true, data, 6) + + +def test_get_layer_heights_agl(): + """Test get_layer_heights with interpolation.""" + heights = np.arange(300, 1200, 100) * units.m + data = heights.m * 0.1 * units.degC + heights, data = get_layer_heights(heights, 500 * units.m, data, with_agl=True) + heights_true = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5]) * units.km + data_true = np.array([30, 40, 50, 60, 70, 80]) * units.degC + assert_array_almost_equal(heights_true, heights, 6) + assert_array_almost_equal(data_true, data, 6) + + +def test_get_layer_heights_agl_bottom_no_interp(): + """Test get_layer_heights with no interpolation and a bottom.""" + heights = np.arange(300, 1200, 100) * units.m + data = heights.m * 0.1 * units.degC + heights, data = get_layer_heights(heights, 500 * units.m, data, with_agl=True, + interpolation=False, bottom=200 * units.m) + heights_true = np.array([0.2, 0.3, 0.4, 0.5, 0.6, 0.7]) * units.km + data_true = np.array([50, 60, 70, 80, 90, 100]) * units.degC + assert_array_almost_equal(heights_true, heights, 6) + assert_array_almost_equal(data_true, data, 6) From 6f08e56049c9be211a42706c41acb0fe23710889 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Tue, 3 Oct 2017 09:09:28 -0600 Subject: [PATCH 4/7] [MNT] Simplify storm relative helicity tests --- metpy/calc/tests/test_kinematics.py | 69 ++++++++++++++--------------- 1 file changed, 33 insertions(+), 36 deletions(-) diff --git a/metpy/calc/tests/test_kinematics.py b/metpy/calc/tests/test_kinematics.py index 23366953a45..cc882a20aa3 100644 --- a/metpy/calc/tests/test_kinematics.py +++ b/metpy/calc/tests/test_kinematics.py @@ -7,8 +7,7 @@ import pytest from metpy.calc import (advection, convergence_vorticity, frontogenesis, geostrophic_wind, - get_wind_components, h_convergence, - montgomery_streamfunction, shearing_deformation, + h_convergence, montgomery_streamfunction, shearing_deformation, shearing_stretching_deformation, storm_relative_helicity, stretching_deformation, total_deformation, v_vorticity) from metpy.constants import g, omega, Re @@ -394,37 +393,35 @@ def test_streamfunc(): assert_almost_equal(msf, 337468.2500 * units('m^2 s^-2'), 4) -def test_helicity(): - """Test function for SRH calculations on an eigth-circle hodograph.""" - pres_test = np.asarray([1013.25, 954.57955706, 898.690770743, - 845.481604002, 794.85264282]) * units('hPa') - hgt_test = hgt_test = np.asarray([0, 500, 1000, 1500, 2000]) - - # Create larger arrays for everything except pressure to make a smoother graph - hgt_int = np.arange(0, 2050, 50) - hgt_int = hgt_int * units('meter') - dir_int = np.arange(180, 272.25, 2.25) - spd_int = np.zeros((hgt_int.shape[0])) - spd_int[:] = 2. - u_int, v_int = get_wind_components(spd_int * units('m/s'), dir_int * units.degree) - - # Interpolate pressure to that graph - pres_int = np.interp(hgt_int, hgt_test, np.log(pres_test.magnitude)) - pres_int = np.exp(pres_int) * units('hPa') - - # Put in the correct value of SRH for a eighth-circle, 2 m/s hodograph - # (SRH = 2 * area under hodo, in this case...) - srh_true_p = (.25 * np.pi * (2 ** 2)) * units('m^2/s^2') - - # Since there's only positive SRH in this case, total SRH will be equal to positive SRH and - # negative SRH will be zero. - srh_true_t = srh_true_p - srh_true_n = 0 * units('m^2/s^2') - p_srh, n_srh, T_srh = storm_relative_helicity(u_int, v_int, pres_int, - hgt_int, 1000 * units('meter'), - bottom=0 * units('meter'), - storm_u=0 * units.knot, - storm_v=0 * units.knot) - assert_almost_equal(p_srh, srh_true_p, 2) - assert_almost_equal(n_srh, srh_true_n, 2) - assert_almost_equal(T_srh, srh_true_t, 2) +def test_storm_relative_helicity_no_storm_motion(): + """Test storm relative helicity with no storm motion and differing input units.""" + u = np.array([0, 20, 10, 0]) * units('m/s') + v = np.array([20, 0, 0, 10]) * units('m/s') + u = u.to('knots') + pressure = np.array([1000, 900, 800, 700]) * units.hPa + heights = np.array([0, 250, 500, 750]) * units.m + + positive_srh, negative_srh, total_srh = storm_relative_helicity(u, v, pressure, heights, + top=750 * units.meters) + + assert_almost_equal(positive_srh, 400. * units('meter ** 2 / second ** 2 '), 6) + assert_almost_equal(negative_srh, -100. * units('meter ** 2 / second ** 2 '), 6) + assert_almost_equal(total_srh, 300. * units('meter ** 2 / second ** 2 '), 6) + + +def test_storm_relative_helicity_storm_motion(): + """Test storm relative helicity with storm motion and differing input units.""" + u = np.array([5, 25, 15, 5]) * units('m/s') + v = np.array([30, 10, 10, 20]) * units('m/s') + u = u.to('knots') + pressure = np.array([1000, 900, 800, 700]) * units.hPa + heights = np.array([0, 250, 500, 750]) * units.m + + pos_srh, neg_srh, total_srh = storm_relative_helicity(u, v, pressure, heights, + top=750 * units.meters, + storm_u=5 * units('m/s'), + storm_v=10 * units('m/s')) + + assert_almost_equal(pos_srh, 400. * units('meter ** 2 / second ** 2 '), 6) + assert_almost_equal(neg_srh, -100. * units('meter ** 2 / second ** 2 '), 6) + assert_almost_equal(total_srh, 300. * units('meter ** 2 / second ** 2 '), 6) From 02d90261df0328322a1c149d6af4dfe4dc2d2eff Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Tue, 3 Oct 2017 13:09:02 -0600 Subject: [PATCH 5/7] Refactor SRH --- metpy/calc/kinematics.py | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/metpy/calc/kinematics.py b/metpy/calc/kinematics.py index 0f68f05b5e4..aa982140121 100644 --- a/metpy/calc/kinematics.py +++ b/metpy/calc/kinematics.py @@ -9,7 +9,7 @@ import numpy as np -from ..calc.tools import get_layer +from ..calc.tools import get_layer_heights from ..cbook import is_string_like, iterable from ..constants import Cp_d, g from ..package_tools import Exporter @@ -539,9 +539,9 @@ def montgomery_streamfunction(height, temperature): @exporter.export -@check_units('[speed]', '[speed]', '[pressure]', '[length]', - '[length]', '[length]', '[speed]', '[speed]') -def storm_relative_helicity(u, v, p, heights, top, bottom=0 * units('meter'), +@check_units('[speed]', '[speed]', '[length]', '[length]', '[length]', + '[speed]', '[speed]') +def storm_relative_helicity(u, v, heights, depth, bottom=0 * units.m, storm_u=0 * units('m/s'), storm_v=0 * units('m/s')): # Partially adapted from similar SharpPy code r"""Calculate storm relative helicity. @@ -561,18 +561,16 @@ def storm_relative_helicity(u, v, p, heights, top, bottom=0 * units('meter'), u component winds v : array-like v component winds - p : array-like - atmospheric pressure heights : array-like - atmospheric heights - top : number - layer top + atmospheric heights, will be converted to AGL + depth : number + depth of the layer bottom : number - layer bottom (default is surface) + height of layer bottom AGL (default is surface) storm_u : number - u component of storm motion + u component of storm motion (default is 0 m/s) storm_v : number - v component of storm motion + v component of storm motion (default is 0 m/s) Returns ------- @@ -580,11 +578,10 @@ def storm_relative_helicity(u, v, p, heights, top, bottom=0 * units('meter'), positive, negative, total storm-relative helicity """ - _, layer_u, layer_v = get_layer(p, u, v, heights=heights, - bottom=bottom, depth=top - bottom) + _, u, v = get_layer_heights(heights, depth, u, v, with_agl=True, bottom=bottom) - storm_relative_u = layer_u - storm_u - storm_relative_v = layer_v - storm_v + storm_relative_u = u - storm_u + storm_relative_v = v - storm_v int_layers = (storm_relative_u[1:] * storm_relative_v[:-1] - storm_relative_u[:-1] * storm_relative_v[1:]) From c585b7a7ce3e3bfed6001de2ad3bf7eeb853ac25 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Tue, 3 Oct 2017 13:09:17 -0600 Subject: [PATCH 6/7] Extend tests for SRH --- metpy/calc/tests/test_kinematics.py | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/metpy/calc/tests/test_kinematics.py b/metpy/calc/tests/test_kinematics.py index cc882a20aa3..bc5ec3f4353 100644 --- a/metpy/calc/tests/test_kinematics.py +++ b/metpy/calc/tests/test_kinematics.py @@ -398,11 +398,10 @@ def test_storm_relative_helicity_no_storm_motion(): u = np.array([0, 20, 10, 0]) * units('m/s') v = np.array([20, 0, 0, 10]) * units('m/s') u = u.to('knots') - pressure = np.array([1000, 900, 800, 700]) * units.hPa heights = np.array([0, 250, 500, 750]) * units.m - positive_srh, negative_srh, total_srh = storm_relative_helicity(u, v, pressure, heights, - top=750 * units.meters) + positive_srh, negative_srh, total_srh = storm_relative_helicity(u, v, heights, + depth=750 * units.meters) assert_almost_equal(positive_srh, 400. * units('meter ** 2 / second ** 2 '), 6) assert_almost_equal(negative_srh, -100. * units('meter ** 2 / second ** 2 '), 6) @@ -414,11 +413,28 @@ def test_storm_relative_helicity_storm_motion(): u = np.array([5, 25, 15, 5]) * units('m/s') v = np.array([30, 10, 10, 20]) * units('m/s') u = u.to('knots') - pressure = np.array([1000, 900, 800, 700]) * units.hPa heights = np.array([0, 250, 500, 750]) * units.m - pos_srh, neg_srh, total_srh = storm_relative_helicity(u, v, pressure, heights, - top=750 * units.meters, + pos_srh, neg_srh, total_srh = storm_relative_helicity(u, v, heights, + depth=750 * units.meters, + storm_u=5 * units('m/s'), + storm_v=10 * units('m/s')) + + assert_almost_equal(pos_srh, 400. * units('meter ** 2 / second ** 2 '), 6) + assert_almost_equal(neg_srh, -100. * units('meter ** 2 / second ** 2 '), 6) + assert_almost_equal(total_srh, 300. * units('meter ** 2 / second ** 2 '), 6) + + +def test_storm_relative_helicity_with_interpolation(): + """Test storm relative helicity with interpolation.""" + u = np.array([-5, 15, 25, 15, -5]) * units('m/s') + v = np.array([40, 20, 10, 10, 30]) * units('m/s') + u = u.to('knots') + heights = np.array([0, 100, 200, 300, 400]) * units.m + + pos_srh, neg_srh, total_srh = storm_relative_helicity(u, v, heights, + bottom=50 * units.meters, + depth=300 * units.meters, storm_u=5 * units('m/s'), storm_v=10 * units('m/s')) From 120eab6c8096554d080a6225aa116b06b6c2f701 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Mon, 9 Oct 2017 13:07:22 -0600 Subject: [PATCH 7/7] Add refactored version of original test. --- metpy/calc/tests/test_kinematics.py | 35 ++++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 3 deletions(-) diff --git a/metpy/calc/tests/test_kinematics.py b/metpy/calc/tests/test_kinematics.py index bc5ec3f4353..2d1c8f86e88 100644 --- a/metpy/calc/tests/test_kinematics.py +++ b/metpy/calc/tests/test_kinematics.py @@ -7,9 +7,10 @@ import pytest from metpy.calc import (advection, convergence_vorticity, frontogenesis, geostrophic_wind, - h_convergence, montgomery_streamfunction, shearing_deformation, - shearing_stretching_deformation, storm_relative_helicity, - stretching_deformation, total_deformation, v_vorticity) + get_wind_components, h_convergence, montgomery_streamfunction, + shearing_deformation, shearing_stretching_deformation, + storm_relative_helicity, stretching_deformation, total_deformation, + v_vorticity) from metpy.constants import g, omega, Re from metpy.testing import assert_almost_equal, assert_array_equal from metpy.units import concatenate, units @@ -441,3 +442,31 @@ def test_storm_relative_helicity_with_interpolation(): assert_almost_equal(pos_srh, 400. * units('meter ** 2 / second ** 2 '), 6) assert_almost_equal(neg_srh, -100. * units('meter ** 2 / second ** 2 '), 6) assert_almost_equal(total_srh, 300. * units('meter ** 2 / second ** 2 '), 6) + + +def test_storm_relative_helicity(): + """Test function for SRH calculations on an eigth-circle hodograph.""" + # Create larger arrays for everything except pressure to make a smoother graph + hgt_int = np.arange(0, 2050, 50) + hgt_int = hgt_int * units('meter') + dir_int = np.arange(180, 272.25, 2.25) + spd_int = np.zeros((hgt_int.shape[0])) + spd_int[:] = 2. + u_int, v_int = get_wind_components(spd_int * units('m/s'), dir_int * units.degree) + + # Put in the correct value of SRH for a eighth-circle, 2 m/s hodograph + # (SRH = 2 * area under hodo, in this case...) + srh_true_p = (.25 * np.pi * (2 ** 2)) * units('m^2/s^2') + + # Since there's only positive SRH in this case, total SRH will be equal to positive SRH and + # negative SRH will be zero. + srh_true_t = srh_true_p + srh_true_n = 0 * units('m^2/s^2') + p_srh, n_srh, T_srh = storm_relative_helicity(u_int, v_int, + hgt_int, 1000 * units('meter'), + bottom=0 * units('meter'), + storm_u=0 * units.knot, + storm_v=0 * units.knot) + assert_almost_equal(p_srh, srh_true_p, 2) + assert_almost_equal(n_srh, srh_true_n, 2) + assert_almost_equal(T_srh, srh_true_t, 2)