forked from Unidata/MetPy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
test_indices.py
217 lines (178 loc) · 9.7 KB
/
test_indices.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
# Copyright (c) 2017,2019 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Test the `indices` module."""
from datetime import datetime
import numpy as np
import xarray as xr
from metpy.calc import (bulk_shear, bunkers_storm_motion, critical_angle,
mean_pressure_weighted, precipitable_water,
significant_tornado, supercell_composite)
from metpy.testing import assert_almost_equal, assert_array_almost_equal, get_upper_air_data
from metpy.units import concatenate, units
def test_precipitable_water():
"""Test precipitable water with observed sounding."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
pw = precipitable_water(data['pressure'], data['dewpoint'], top=400 * units.hPa)
truth = 22.60430651 * units.millimeters
assert_array_almost_equal(pw, truth, 4)
def test_precipitable_water_no_bounds():
"""Test precipitable water with observed sounding and no bounds given."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
dewpoint = data['dewpoint']
pressure = data['pressure']
inds = pressure >= 400 * units.hPa
pw = precipitable_water(pressure[inds], dewpoint[inds])
truth = 22.60430651 * units.millimeters
assert_array_almost_equal(pw, truth, 4)
def test_precipitable_water_bound_error():
"""Test with no top bound given and data that produced floating point issue #596."""
pressure = np.array([993., 978., 960.5, 927.6, 925., 895.8, 892., 876., 45.9, 39.9, 36.,
36., 34.3]) * units.hPa
dewpoint = np.array([25.5, 24.1, 23.1, 21.2, 21.1, 19.4, 19.2, 19.2, -87.1, -86.5, -86.5,
-86.5, -88.1]) * units.degC
pw = precipitable_water(pressure, dewpoint)
truth = 89.86846252697836 * units('millimeters')
assert_almost_equal(pw, truth, 5)
def test_precipitable_water_nans():
"""Test that PW returns appropriate number if NaNs are present."""
pressure = np.array([1001, 1000, 997, 977.9, 977, 957, 937.8, 925, 906, 899.3, 887, 862.5,
854, 850, 800, 793.9, 785, 777, 771, 762, 731.8, 726, 703, 700, 655,
630, 621.2, 602, 570.7, 548, 546.8, 539, 513, 511, 485, 481, 468,
448, 439, 424, 420, 412]) * units.hPa
dewpoint = np.array([-25.1, -26.1, -26.8, np.nan, -27.3, -28.2, np.nan, -27.2, -26.6,
np.nan, -27.4, np.nan, -23.5, -23.5, -25.1, np.nan, -22.9, -17.8,
-16.6, np.nan, np.nan, -16.4, np.nan, -18.5, -21., -23.7, np.nan,
-28.3, np.nan, -32.6, np.nan, -33.8, -35., -35.1, -38.1, -40.,
-43.3, -44.6, -46.4, -47., -49.2, -50.7]) * units.degC
pw = precipitable_water(pressure, dewpoint)
truth = 4.003660322395436 * units.mm
assert_almost_equal(pw, truth, 5)
def test_mean_pressure_weighted():
"""Test pressure-weighted mean wind function with vertical interpolation."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
u, v = mean_pressure_weighted(data['pressure'],
data['u_wind'],
data['v_wind'],
height=data['height'],
depth=6000 * units('meter'))
assert_almost_equal(u, 6.0208700094534775 * units('m/s'), 7)
assert_almost_equal(v, 7.966031839967931 * units('m/s'), 7)
def test_mean_pressure_weighted_temperature():
"""Test pressure-weighted mean temperature function with vertical interpolation."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
T, = mean_pressure_weighted(data['pressure'],
data['temperature'],
height=data['height'],
depth=6000 * units('meter'))
assert_almost_equal(T, 281.535035296836 * units('kelvin'), 7)
def test_mean_pressure_weighted_elevated():
"""Test pressure-weighted mean wind function with a base above the surface."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
u, v = mean_pressure_weighted(data['pressure'],
data['u_wind'],
data['v_wind'],
height=data['height'],
depth=3000 * units('meter'),
bottom=data['height'][0] + 3000 * units('meter'))
assert_almost_equal(u, 8.270829843626476 * units('m/s'), 7)
assert_almost_equal(v, 1.7392601775853547 * units('m/s'), 7)
def test_precipitable_water_xarray():
"""Test precipitable water with xarray input."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
press = xr.DataArray(data['pressure'], attrs={'units': str(data['pressure'].units)})
dewp = xr.DataArray(data['dewpoint'], dims=('press',), coords=(press,))
pw = precipitable_water(press, dewp, top=400 * units.hPa)
truth = 22.60430651 * units.millimeters
assert_almost_equal(pw, truth)
def test_bunkers_motion():
"""Test Bunkers storm motion with observed sounding."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
motion = concatenate(bunkers_storm_motion(data['pressure'],
data['u_wind'], data['v_wind'],
data['height']))
truth = [1.4537892577864744, 2.0169333025630616, 10.587950761120482, 13.915130377372801,
6.0208700094534775, 7.9660318399679308] * units('m/s')
assert_almost_equal(motion.flatten(), truth, 8)
def test_bulk_shear():
"""Test bulk shear with observed sounding."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
u, v = bulk_shear(data['pressure'], data['u_wind'],
data['v_wind'], height=data['height'],
depth=6000 * units('meter'))
truth = [29.899581266946115, -14.389225800205509] * units('knots')
assert_almost_equal(u.to('knots'), truth[0], 8)
assert_almost_equal(v.to('knots'), truth[1], 8)
def test_bulk_shear_no_depth():
"""Test bulk shear with observed sounding and no depth given. Issue #568."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
u, v = bulk_shear(data['pressure'], data['u_wind'],
data['v_wind'], height=data['height'])
truth = [20.225018939, 22.602359692] * units('knots')
assert_almost_equal(u.to('knots'), truth[0], 8)
assert_almost_equal(v.to('knots'), truth[1], 8)
def test_bulk_shear_elevated():
"""Test bulk shear with observed sounding and a base above the surface."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
u, v = bulk_shear(data['pressure'], data['u_wind'],
data['v_wind'], height=data['height'],
bottom=data['height'][0] + 3000 * units('meter'),
depth=3000 * units('meter'))
truth = [0.9655943923302139, -3.8405428777944466] * units('m/s')
assert_almost_equal(u, truth[0], 8)
assert_almost_equal(v, truth[1], 8)
def test_supercell_composite():
"""Test supercell composite function."""
mucape = [2000., 1000., 500., 2000.] * units('J/kg')
esrh = [400., 150., 45., 45.] * units('m^2/s^2')
ebwd = [30., 15., 5., 5.] * units('m/s')
truth = [16., 2.25, 0., 0.]
supercell_comp = supercell_composite(mucape, esrh, ebwd)
assert_array_almost_equal(supercell_comp, truth, 5)
def test_supercell_composite_scalar():
"""Test supercell composite function with a single value."""
mucape = 2000. * units('J/kg')
esrh = 400. * units('m^2/s^2')
ebwd = 30. * units('m/s')
truth = 16.
supercell_comp = supercell_composite(mucape, esrh, ebwd)
assert_almost_equal(supercell_comp, truth, 6)
def test_sigtor():
"""Test significant tornado parameter function."""
sbcape = [2000., 2000., 2000., 2000., 3000, 4000] * units('J/kg')
sblcl = [3000., 1500., 500., 1500., 1500, 800] * units('meter')
srh1 = [200., 200., 200., 200., 300, 400] * units('m^2/s^2')
shr6 = [20., 5., 20., 35., 20., 35] * units('m/s')
truth = [0., 0, 1.777778, 1.333333, 2., 10.666667]
sigtor = significant_tornado(sbcape, sblcl, srh1, shr6)
assert_almost_equal(sigtor, truth, 6)
def test_sigtor_scalar():
"""Test significant tornado parameter function with a single value."""
sbcape = 4000 * units('J/kg')
sblcl = 800 * units('meter')
srh1 = 400 * units('m^2/s^2')
shr6 = 35 * units('m/s')
truth = 10.666667
sigtor = significant_tornado(sbcape, sblcl, srh1, shr6)
assert_almost_equal(sigtor, truth, 6)
def test_critical_angle():
"""Test critical angle with observed sounding."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
ca = critical_angle(data['pressure'], data['u_wind'],
data['v_wind'], data['height'],
u_storm=0 * units('m/s'), v_storm=0 * units('m/s'))
truth = [140.0626637513269] * units('degrees')
assert_almost_equal(ca, truth, 8)
def test_critical_angle_units():
"""Test critical angle with observed sounding and different storm motion units."""
data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')
# Set storm motion in m/s
ca_ms = critical_angle(data['pressure'], data['u_wind'],
data['v_wind'], data['height'],
u_storm=10 * units('m/s'), v_storm=10 * units('m/s'))
# Set same storm motion in kt and m/s
ca_kt_ms = critical_angle(data['pressure'], data['u_wind'],
data['v_wind'], data['height'],
u_storm=10 * units('m/s'), v_storm=19.4384449244 * units('kt'))
# Make sure the resulting critical angles are equal
assert_almost_equal(ca_ms, ca_kt_ms, 8)