|
1
|
|
|
# Copyright (c) 2008-2015 MetPy Developers. |
|
2
|
|
|
# Distributed under the terms of the BSD 3-Clause License. |
|
3
|
|
|
# SPDX-License-Identifier: BSD-3-Clause |
|
4
|
|
|
"""Test the `thermo` module.""" |
|
5
|
|
|
|
|
6
|
|
|
import numpy as np |
|
7
|
|
|
import pytest |
|
8
|
|
|
|
|
9
|
|
|
from metpy.calc import (cape_cin, density, dewpoint, dewpoint_rh, dry_lapse, el, |
|
10
|
|
|
equivalent_potential_temperature, |
|
11
|
|
|
isentropic_interpolation, |
|
12
|
|
|
lcl, lfc, mixing_ratio, |
|
13
|
|
|
mixing_ratio_from_specific_humidity, moist_lapse, |
|
14
|
|
|
most_unstable_parcel, parcel_profile, potential_temperature, |
|
15
|
|
|
psychrometric_vapor_pressure_wet, |
|
16
|
|
|
relative_humidity_from_mixing_ratio, |
|
17
|
|
|
relative_humidity_from_specific_humidity, |
|
18
|
|
|
relative_humidity_wet_psychrometric, |
|
19
|
|
|
saturation_mixing_ratio, |
|
20
|
|
|
saturation_vapor_pressure, vapor_pressure, |
|
21
|
|
|
virtual_potential_temperature, virtual_temperature) |
|
22
|
|
|
|
|
23
|
|
|
from metpy.calc.thermo import _find_append_zero_crossings |
|
24
|
|
|
from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan |
|
25
|
|
|
from metpy.units import units |
|
26
|
|
|
|
|
27
|
|
|
|
|
28
|
|
|
def test_potential_temperature(): |
|
29
|
|
|
"""Test potential_temperature calculation.""" |
|
30
|
|
|
temp = np.array([278., 283., 291., 298.]) * units.kelvin |
|
31
|
|
|
pres = np.array([900., 500., 300., 100.]) * units.mbar |
|
32
|
|
|
real_th = np.array([286.493, 344.961, 410.4335, 575.236]) * units.kelvin |
|
33
|
|
|
assert_array_almost_equal(potential_temperature(pres, temp), real_th, 3) |
|
34
|
|
|
|
|
35
|
|
|
|
|
36
|
|
|
def test_scalar(): |
|
37
|
|
|
"""Test potential_temperature accepts scalar values.""" |
|
38
|
|
|
assert_almost_equal(potential_temperature(1000. * units.mbar, 293. * units.kelvin), |
|
39
|
|
|
293. * units.kelvin, 4) |
|
40
|
|
|
assert_almost_equal(potential_temperature(800. * units.mbar, 293. * units.kelvin), |
|
41
|
|
|
312.2828 * units.kelvin, 4) |
|
42
|
|
|
|
|
43
|
|
|
|
|
44
|
|
|
def test_fahrenheit(): |
|
45
|
|
|
"""Test that potential_temperature handles temperature values in Fahrenheit.""" |
|
46
|
|
|
assert_almost_equal(potential_temperature(800. * units.mbar, 68. * units.degF), |
|
47
|
|
|
(312.444 * units.kelvin).to(units.degF), 2) |
|
48
|
|
|
|
|
49
|
|
|
|
|
50
|
|
|
def test_pot_temp_inhg(): |
|
51
|
|
|
"""Test that potential_temperature can handle pressure not in mb (issue #165).""" |
|
52
|
|
|
assert_almost_equal(potential_temperature(29.92 * units.inHg, 29 * units.degC), |
|
53
|
|
|
301.019735 * units.kelvin, 4) |
|
54
|
|
|
|
|
55
|
|
|
|
|
56
|
|
|
def test_dry_lapse(): |
|
57
|
|
|
"""Test dry_lapse calculation.""" |
|
58
|
|
|
levels = np.array([1000, 900, 864.89]) * units.mbar |
|
59
|
|
|
temps = dry_lapse(levels, 303.15 * units.kelvin) |
|
60
|
|
|
assert_array_almost_equal(temps, |
|
61
|
|
|
np.array([303.15, 294.16, 290.83]) * units.kelvin, 2) |
|
62
|
|
|
|
|
63
|
|
|
|
|
64
|
|
|
def test_dry_lapse_2_levels(): |
|
65
|
|
|
"""Test dry_lapse calculation when given only two levels.""" |
|
66
|
|
|
temps = dry_lapse(np.array([1000., 500.]) * units.mbar, 293. * units.kelvin) |
|
67
|
|
|
assert_array_almost_equal(temps, [293., 240.3723] * units.kelvin, 4) |
|
68
|
|
|
|
|
69
|
|
|
|
|
70
|
|
|
def test_moist_lapse(): |
|
71
|
|
|
"""Test moist_lapse calculation.""" |
|
72
|
|
|
temp = moist_lapse(np.array([1000., 800., 600., 500., 400.]) * units.mbar, |
|
73
|
|
|
293. * units.kelvin) |
|
74
|
|
|
true_temp = np.array([293, 284.64, 272.81, 264.42, 252.91]) * units.kelvin |
|
75
|
|
|
assert_array_almost_equal(temp, true_temp, 2) |
|
76
|
|
|
|
|
77
|
|
|
|
|
78
|
|
|
def test_moist_lapse_degc(): |
|
79
|
|
|
"""Test moist_lapse with Celsius temperatures.""" |
|
80
|
|
|
temp = moist_lapse(np.array([1000., 800., 600., 500., 400.]) * units.mbar, |
|
81
|
|
|
19.85 * units.degC) |
|
82
|
|
|
true_temp = np.array([293, 284.64, 272.81, 264.42, 252.91]) * units.kelvin |
|
83
|
|
|
assert_array_almost_equal(temp, true_temp, 2) |
|
84
|
|
|
|
|
85
|
|
|
|
|
86
|
|
|
def test_parcel_profile(): |
|
87
|
|
|
"""Test parcel profile calculation.""" |
|
88
|
|
|
levels = np.array([1000., 900., 800., 700., 600., 500., 400.]) * units.mbar |
|
89
|
|
|
true_prof = np.array([303.15, 294.16, 288.026, 283.073, 277.058, 269.402, |
|
90
|
|
|
258.966]) * units.kelvin |
|
91
|
|
|
|
|
92
|
|
|
prof = parcel_profile(levels, 30. * units.degC, 20. * units.degC) |
|
93
|
|
|
assert_array_almost_equal(prof, true_prof, 2) |
|
94
|
|
|
|
|
95
|
|
|
|
|
96
|
|
|
def test_parcel_profile_saturated(): |
|
97
|
|
|
"""Test parcel_profile works when LCL in levels (issue #232).""" |
|
98
|
|
|
levels = np.array([1000., 700., 500.]) * units.mbar |
|
99
|
|
|
true_prof = np.array([296.95, 284.381, 271.123]) * units.kelvin |
|
100
|
|
|
|
|
101
|
|
|
prof = parcel_profile(levels, 23.8 * units.degC, 23.8 * units.degC) |
|
102
|
|
|
assert_array_almost_equal(prof, true_prof, 2) |
|
103
|
|
|
|
|
104
|
|
|
|
|
105
|
|
|
def test_sat_vapor_pressure(): |
|
106
|
|
|
"""Test saturation_vapor_pressure calculation.""" |
|
107
|
|
|
temp = np.array([5., 10., 18., 25.]) * units.degC |
|
108
|
|
|
real_es = np.array([8.72, 12.27, 20.63, 31.67]) * units.mbar |
|
109
|
|
|
assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 2) |
|
110
|
|
|
|
|
111
|
|
|
|
|
112
|
|
|
def test_sat_vapor_pressure_scalar(): |
|
113
|
|
|
"""Test saturation_vapor_pressure handles scalar values.""" |
|
114
|
|
|
es = saturation_vapor_pressure(0 * units.degC) |
|
115
|
|
|
assert_almost_equal(es, 6.112 * units.mbar, 3) |
|
116
|
|
|
|
|
117
|
|
|
|
|
118
|
|
|
def test_sat_vapor_pressure_fahrenheit(): |
|
119
|
|
|
"""Test saturation_vapor_pressure handles temperature in Fahrenheit.""" |
|
120
|
|
|
temp = np.array([50., 68.]) * units.degF |
|
121
|
|
|
real_es = np.array([12.2717, 23.3695]) * units.mbar |
|
122
|
|
|
assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 4) |
|
123
|
|
|
|
|
124
|
|
|
|
|
125
|
|
|
def test_basic_dewpoint_rh(): |
|
126
|
|
|
"""Test dewpoint_rh function.""" |
|
127
|
|
|
temp = np.array([30., 25., 10., 20., 25.]) * units.degC |
|
128
|
|
|
rh = np.array([30., 45., 55., 80., 85.]) / 100. |
|
129
|
|
|
|
|
130
|
|
|
real_td = np.array([11, 12, 1, 16, 22]) * units.degC |
|
131
|
|
|
assert_array_almost_equal(real_td, dewpoint_rh(temp, rh), 0) |
|
132
|
|
|
|
|
133
|
|
|
|
|
134
|
|
|
def test_scalar_dewpoint_rh(): |
|
135
|
|
|
"""Test dewpoint_rh with scalar values.""" |
|
136
|
|
|
td = dewpoint_rh(10.6 * units.degC, 0.37) |
|
137
|
|
|
assert_almost_equal(td, 26. * units.degF, 0) |
|
138
|
|
|
|
|
139
|
|
|
|
|
140
|
|
|
def test_dewpoint(): |
|
141
|
|
|
"""Test dewpoint calculation.""" |
|
142
|
|
|
assert_almost_equal(dewpoint(6.112 * units.mbar), 0. * units.degC, 2) |
|
143
|
|
|
|
|
144
|
|
|
|
|
145
|
|
|
def test_dewpoint_weird_units(): |
|
146
|
|
|
"""Test dewpoint using non-standard units. |
|
147
|
|
|
|
|
148
|
|
|
Revealed from odd dimensionless units and ending up using numpy.ma math |
|
149
|
|
|
functions instead of numpy ones. |
|
150
|
|
|
""" |
|
151
|
|
|
assert_almost_equal(dewpoint(15825.6 * units('g * mbar / kg')), |
|
152
|
|
|
13.8564 * units.degC, 4) |
|
153
|
|
|
|
|
154
|
|
|
|
|
155
|
|
|
def test_mixing_ratio(): |
|
156
|
|
|
"""Test mixing ratio calculation.""" |
|
157
|
|
|
p = 998. * units.mbar |
|
158
|
|
|
e = 73.75 * units.mbar |
|
159
|
|
|
assert_almost_equal(mixing_ratio(e, p), 0.04963, 2) |
|
160
|
|
|
|
|
161
|
|
|
|
|
162
|
|
|
def test_vapor_pressure(): |
|
163
|
|
|
"""Test vapor pressure calculation.""" |
|
164
|
|
|
assert_almost_equal(vapor_pressure(998. * units.mbar, 0.04963), |
|
165
|
|
|
73.74925 * units.mbar, 5) |
|
166
|
|
|
|
|
167
|
|
|
|
|
168
|
|
|
def test_lcl(): |
|
169
|
|
|
"""Test LCL calculation.""" |
|
170
|
|
|
lcl_pressure, lcl_temperature = lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC) |
|
171
|
|
|
assert_almost_equal(lcl_pressure, 864.761 * units.mbar, 2) |
|
172
|
|
|
assert_almost_equal(lcl_temperature, 17.676 * units.degC, 2) |
|
173
|
|
|
|
|
174
|
|
|
|
|
175
|
|
|
def test_lcl_convergence(): |
|
176
|
|
|
"""Test LCL calculation convergence failure.""" |
|
177
|
|
|
with pytest.raises(RuntimeError): |
|
178
|
|
|
lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC, max_iters=2) |
|
179
|
|
View Code Duplication |
|
|
|
|
|
|
|
180
|
|
|
|
|
181
|
|
|
def test_lfc_basic(): |
|
182
|
|
|
"""Test LFC calculation.""" |
|
183
|
|
|
levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar |
|
184
|
|
|
temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -49.]) * units.celsius |
|
185
|
|
|
dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius |
|
186
|
|
|
l = lfc(levels, temperatures, dewpoints) |
|
187
|
|
|
assert_almost_equal(l[0], 727.468 * units.mbar, 2) |
|
188
|
|
|
assert_almost_equal(l[1], 9.705 * units.celsius, 2) |
|
189
|
|
|
|
|
190
|
|
|
|
|
191
|
|
|
def test_no_lfc(): |
|
192
|
|
|
"""Test LFC calculation when there is no LFC in the data.""" |
|
193
|
|
|
levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar |
|
194
|
|
|
temperatures = np.array([22.2, 17.4, 14.6, 1.4, -17.6, -39.4, -52.5]) * units.celsius |
|
195
|
|
|
dewpoints = np.array([9., 4.3, -21.2, -26.7, -31., -53.3, -66.7]) * units.celsius |
|
196
|
|
|
lfc_pressure, lfc_temperature = lfc(levels, temperatures, dewpoints) |
|
197
|
|
|
assert assert_nan(lfc_pressure, levels.units) |
|
198
|
|
|
assert assert_nan(lfc_temperature, temperatures.units) |
|
199
|
|
|
|
|
200
|
|
|
|
|
201
|
|
|
def test_lfc_inversion(): |
|
202
|
|
|
"""Test LFC when there is an inversion to be sure we don't pick that.""" |
|
203
|
|
|
levels = np.array([963., 789., 782.3, 754.8, 728.1, 727., 700., |
|
204
|
|
|
571., 450., 300., 248.]) * units.mbar |
|
205
|
|
|
temperatures = np.array([25.4, 18.4, 17.8, 15.4, 12.9, 12.8, |
|
206
|
|
|
10., -3.9, -16.3, -41.1, -51.5]) * units.celsius |
|
207
|
|
|
dewpoints = np.array([20.4, 0.4, -0.5, -4.3, -8., -8.2, -9., |
|
208
|
|
|
-23.9, -33.3, -54.1, -63.5]) * units.celsius |
|
209
|
|
|
l = lfc(levels, temperatures, dewpoints) |
|
210
|
|
|
assert_almost_equal(l[0], 706.0103 * units.mbar, 2) |
|
211
|
|
|
assert_almost_equal(l[1], 10.6232 * units.celsius, 2) |
|
212
|
|
|
|
|
213
|
|
|
|
|
214
|
|
|
def test_lfc_equals_lcl(): |
|
215
|
|
|
"""Test LFC when there is no cap and the lfc is equal to the lcl.""" |
|
216
|
|
|
levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1, |
|
217
|
|
|
748., 732.2, 700., 654.8]) * units.mbar |
|
218
|
|
|
temperatures = np.array([29.4, 28.7, 25.2, 22.4, 19.4, 16.8, |
|
219
|
|
|
14.3, 13.2, 12.6, 11.4, 7.1]) * units.celsius |
|
220
|
|
|
dewpoints = np.array([18.4, 18.1, 16.6, 15.4, 13.2, 11.4, 9.6, |
|
221
|
|
|
8.8, 0., -18.6, -22.9]) * units.celsius |
|
222
|
|
|
l = lfc(levels, temperatures, dewpoints) |
|
223
|
|
|
assert_almost_equal(l[0], 777.0333 * units.mbar, 2) |
|
224
|
|
|
assert_almost_equal(l[1], 15.8714 * units.celsius, 2) |
|
225
|
|
|
|
|
226
|
|
|
|
|
227
|
|
|
def test_lfc_sfc_precision(): |
|
228
|
|
|
"""Test LFC when there are precision issues with the parcel path.""" |
|
229
|
|
|
levels = np.array([839., 819.4, 816., 807., 790.7, 763., 736.2, |
|
230
|
|
|
722., 710.1, 700.]) * units.mbar |
|
231
|
|
|
temperatures = np.array([20.6, 22.3, 22.6, 22.2, 20.9, 18.7, 16.4, |
|
232
|
|
|
15.2, 13.9, 12.8]) * units.celsius |
|
233
|
|
|
dewpoints = np.array([10.6, 8., 7.6, 6.2, 5.7, 4.7, 3.7, 3.2, 3., 2.8]) * units.celsius |
|
234
|
|
|
l = lfc(levels, temperatures, dewpoints) |
|
235
|
|
|
assert assert_nan(l[0], levels.units) |
|
236
|
|
|
assert assert_nan(l[1], temperatures.units) |
|
237
|
|
|
|
|
238
|
|
|
|
|
239
|
|
|
def test_saturation_mixing_ratio(): |
|
240
|
|
|
"""Test saturation mixing ratio calculation.""" |
|
241
|
|
|
p = 999. * units.mbar |
|
242
|
|
|
t = 288. * units.kelvin |
|
243
|
|
|
assert_almost_equal(saturation_mixing_ratio(p, t), .01068, 3) |
|
244
|
|
|
|
|
245
|
|
|
|
|
246
|
|
|
def test_equivalent_potential_temperature(): |
|
247
|
|
|
"""Test equivalent potential temperature calculation.""" |
|
248
|
|
|
p = 999. * units.mbar |
|
249
|
|
|
t = 288. * units.kelvin |
|
250
|
|
|
ept = equivalent_potential_temperature(p, t) |
|
251
|
|
|
assert_almost_equal(ept, 315.9548 * units.kelvin, 3) |
|
252
|
|
|
|
|
253
|
|
|
|
|
254
|
|
|
def test_virtual_temperature(): |
|
255
|
|
|
"""Test virtual temperature calculation.""" |
|
256
|
|
|
t = 288. * units.kelvin |
|
257
|
|
|
qv = .0016 # kg/kg |
|
258
|
|
|
tv = virtual_temperature(t, qv) |
|
259
|
|
|
assert_almost_equal(tv, 288.2796 * units.kelvin, 3) |
|
260
|
|
|
|
|
261
|
|
|
|
|
262
|
|
|
def test_virtual_potential_temperature(): |
|
263
|
|
|
"""Test virtual potential temperature calculation.""" |
|
264
|
|
|
p = 999. * units.mbar |
|
265
|
|
|
t = 288. * units.kelvin |
|
266
|
|
|
qv = .0016 # kg/kg |
|
267
|
|
|
theta_v = virtual_potential_temperature(p, t, qv) |
|
268
|
|
|
assert_almost_equal(theta_v, 288.3620 * units.kelvin, 3) |
|
269
|
|
|
|
|
270
|
|
|
|
|
271
|
|
|
def test_density(): |
|
272
|
|
|
"""Test density calculation.""" |
|
273
|
|
|
p = 999. * units.mbar |
|
274
|
|
|
t = 288. * units.kelvin |
|
275
|
|
|
qv = .0016 # kg/kg |
|
276
|
|
|
rho = density(p, t, qv).to(units.kilogram / units.meter ** 3) |
|
277
|
|
|
assert_almost_equal(rho, 1.2072 * (units.kilogram / units.meter ** 3), 3) |
|
278
|
|
|
|
|
279
|
|
|
|
|
280
|
|
|
def test_el(): |
|
281
|
|
|
"""Test equilibrium layer calculation.""" |
|
282
|
|
|
levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar |
|
283
|
|
|
temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius |
|
284
|
|
|
dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius |
|
285
|
|
|
el_pressure, el_temperature = el(levels, temperatures, dewpoints) |
|
286
|
|
|
assert_almost_equal(el_pressure, 520.8700 * units.mbar, 3) |
|
287
|
|
|
assert_almost_equal(el_temperature, -11.7027 * units.degC, 3) |
|
288
|
|
|
|
|
289
|
|
|
|
|
290
|
|
|
def test_no_el(): |
|
291
|
|
|
"""Test equilibrium layer calculation when there is no EL in the data.""" |
|
292
|
|
|
levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar |
|
293
|
|
|
temperatures = np.array([22.2, 17.4, 14.6, 1.4, -17.6, -39.4, -52.5]) * units.celsius |
|
294
|
|
|
dewpoints = np.array([19., 14.3, -11.2, -16.7, -21., -43.3, -56.7]) * units.celsius |
|
295
|
|
|
el_pressure, el_temperature = el(levels, temperatures, dewpoints) |
|
296
|
|
|
assert assert_nan(el_pressure, levels.units) |
|
297
|
|
|
assert assert_nan(el_temperature, temperatures.units) |
|
298
|
|
|
|
|
299
|
|
|
|
|
300
|
|
|
def test_no_el_multi_crossing(): |
|
301
|
|
|
"""Test el calculation with no el and severel parcel path-profile crossings.""" |
|
302
|
|
|
levels = np.array([918., 911., 880., 873.9, 850., 848., 843.5, 818., 813.8, 785., |
|
303
|
|
|
773., 763., 757.5, 730.5, 700., 679., 654.4, 645., |
|
304
|
|
|
643.9]) * units.mbar |
|
305
|
|
|
temperatures = np.array([24.2, 22.8, 19.6, 19.1, 17., 16.8, 16.5, 15., 14.9, 14.4, 16.4, |
|
306
|
|
|
16.2, 15.7, 13.4, 10.6, 8.4, 5.7, 4.6, 4.5]) * units.celsius |
|
307
|
|
|
dewpoints = np.array([19.5, 17.8, 16.7, 16.5, 15.8, 15.7, 15.3, 13.1, 12.9, 11.9, 6.4, |
|
308
|
|
|
3.2, 2.6, -0.6, -4.4, -6.6, -9.3, -10.4, -10.5]) * units.celsius |
|
309
|
|
|
el_pressure, el_temperature = el(levels, temperatures, dewpoints) |
|
310
|
|
|
assert assert_nan(el_pressure, levels.units) |
|
311
|
|
|
assert assert_nan(el_temperature, temperatures.units) |
|
312
|
|
|
|
|
313
|
|
|
|
|
314
|
|
|
def test_el_lfc_equals_lcl(): |
|
315
|
|
|
"""Test equilibrium layer calculation when the lfc equals the lcl.""" |
|
316
|
|
|
levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1, 748., |
|
317
|
|
|
732.3, 700., 654.8, 606.8, 562.4, 501.8, 500., 482., |
|
318
|
|
|
400., 393.3, 317.1, 307., 300., 252.7, 250., 200., |
|
319
|
|
|
199.3, 197., 190., 172., 156.6, 150., 122.9, 112., |
|
320
|
|
|
106.2, 100.]) * units.mbar |
|
321
|
|
|
temperatures = np.array([29.4, 28.7, 25.2, 22.4, 19.4, 16.8, 14.3, |
|
322
|
|
|
13.2, 12.6, 11.4, 7.1, 2.2, -2.7, -10.1, |
|
323
|
|
|
-10.3, -12.4, -23.3, -24.4, -38., -40.1, -41.1, |
|
324
|
|
|
-49.8, -50.3, -59.1, -59.1, -59.3, -59.7, -56.3, |
|
325
|
|
|
-56.9, -57.1, -59.1, -60.1, -58.6, -56.9]) * units.celsius |
|
326
|
|
|
dewpoints = np.array([18.4, 18.1, 16.6, 15.4, 13.2, 11.4, 9.6, 8.8, 0., |
|
327
|
|
|
-18.6, -22.9, -27.8, -32.7, -40.1, -40.3, -42.4, -53.3, |
|
328
|
|
|
-54.4, -68., -70.1, -70., -70., -70., -70., -70., -70., |
|
329
|
|
|
-70., -70., -70., -70., -70., -70., -70., -70.]) * units.celsius |
|
330
|
|
|
el_pressure, el_temperature = el(levels, temperatures, dewpoints) |
|
331
|
|
|
assert_almost_equal(el_pressure, 175.8684 * units.mbar, 3) |
|
332
|
|
|
assert_almost_equal(el_temperature, -57.0307 * units.degC, 3) |
|
333
|
|
|
|
|
334
|
|
|
|
|
335
|
|
|
def test_wet_psychrometric_vapor_pressure(): |
|
336
|
|
|
"""Test calculation of vapor pressure from wet and dry bulb temperatures.""" |
|
337
|
|
|
p = 1013.25 * units.mbar |
|
338
|
|
|
dry_bulb_temperature = 20. * units.degC |
|
339
|
|
|
wet_bulb_temperature = 18. * units.degC |
|
340
|
|
|
psychrometric_vapor_pressure = psychrometric_vapor_pressure_wet(dry_bulb_temperature, |
|
341
|
|
|
wet_bulb_temperature, p) |
|
342
|
|
|
assert_almost_equal(psychrometric_vapor_pressure, 19.3673 * units.mbar, 3) |
|
343
|
|
|
|
|
344
|
|
|
|
|
345
|
|
|
def test_wet_psychrometric_rh(): |
|
346
|
|
|
"""Test calculation of relative humidity from wet and dry bulb temperatures.""" |
|
347
|
|
|
p = 1013.25 * units.mbar |
|
348
|
|
|
dry_bulb_temperature = 20. * units.degC |
|
349
|
|
|
wet_bulb_temperature = 18. * units.degC |
|
350
|
|
|
psychrometric_rh = relative_humidity_wet_psychrometric(dry_bulb_temperature, |
|
351
|
|
|
wet_bulb_temperature, p) |
|
352
|
|
|
assert_almost_equal(psychrometric_rh, 82.8747 * units.percent, 3) |
|
353
|
|
|
|
|
354
|
|
|
|
|
355
|
|
|
def test_wet_psychrometric_rh_kwargs(): |
|
356
|
|
|
"""Test calculation of relative humidity from wet and dry bulb temperatures.""" |
|
357
|
|
|
p = 1013.25 * units.mbar |
|
358
|
|
|
dry_bulb_temperature = 20. * units.degC |
|
359
|
|
|
wet_bulb_temperature = 18. * units.degC |
|
360
|
|
|
coeff = 6.1e-4 / units.kelvin |
|
361
|
|
|
psychrometric_rh = relative_humidity_wet_psychrometric(dry_bulb_temperature, |
|
362
|
|
|
wet_bulb_temperature, p, |
|
363
|
|
|
psychrometer_coefficient=coeff) |
|
364
|
|
View Code Duplication |
assert_almost_equal(psychrometric_rh, 82.9701 * units.percent, 3) |
|
|
|
|
|
|
365
|
|
|
|
|
366
|
|
|
|
|
367
|
|
|
def test_rh_mixing_ratio(): |
|
368
|
|
|
"""Tests relative humidity from mixing ratio.""" |
|
369
|
|
|
p = 1013.25 * units.mbar |
|
370
|
|
|
temperature = 20. * units.degC |
|
371
|
|
|
w = 0.012 |
|
372
|
|
|
rh = relative_humidity_from_mixing_ratio(w, temperature, p) |
|
373
|
|
|
assert_almost_equal(rh, 81.7219 * units.percent, 3) |
|
374
|
|
|
|
|
375
|
|
|
|
|
376
|
|
|
def test_mixing_ratio_from_specific_humidity(): |
|
377
|
|
|
"""Tests mixing ratio from specific humidity.""" |
|
378
|
|
|
q = 0.012 |
|
379
|
|
|
w = mixing_ratio_from_specific_humidity(q) |
|
380
|
|
|
assert_almost_equal(w, 0.01215, 3) |
|
381
|
|
|
|
|
382
|
|
|
|
|
383
|
|
|
def test_rh_specific_humidity(): |
|
384
|
|
|
"""Tests relative humidity from specific humidity.""" |
|
385
|
|
|
p = 1013.25 * units.mbar |
|
386
|
|
View Code Duplication |
temperature = 20. * units.degC |
|
|
|
|
|
|
387
|
|
|
q = 0.012 |
|
388
|
|
|
rh = relative_humidity_from_specific_humidity(q, temperature, p) |
|
389
|
|
|
assert_almost_equal(rh, 82.7145 * units.percent, 3) |
|
390
|
|
|
|
|
391
|
|
|
|
|
392
|
|
|
def test_cape_cin(): |
|
393
|
|
|
"""Tests the basic CAPE and CIN calculation.""" |
|
394
|
|
|
p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar |
|
395
|
|
|
temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius |
|
396
|
|
|
dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius |
|
397
|
|
|
parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) |
|
398
|
|
|
cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) |
|
399
|
|
|
assert_almost_equal(cape, 58.0368212 * units('joule / kilogram'), 6) |
|
400
|
|
|
assert_almost_equal(cin, -89.8073512 * units('joule / kilogram'), 6) |
|
401
|
|
|
|
|
402
|
|
|
|
|
403
|
|
|
def test_cape_cin_no_el(): |
|
404
|
|
|
"""Tests that CAPE works with no EL.""" |
|
405
|
|
|
p = np.array([959., 779.2, 751.3, 724.3]) * units.mbar |
|
406
|
|
|
temperature = np.array([22.2, 14.6, 12., 9.4]) * units.celsius |
|
407
|
|
|
dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius |
|
408
|
|
|
parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC') |
|
409
|
|
|
cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) |
|
410
|
|
View Code Duplication |
assert_almost_equal(cape, 0.08750805 * units('joule / kilogram'), 6) |
|
|
|
|
|
|
411
|
|
|
assert_almost_equal(cin, -89.8073512 * units('joule / kilogram'), 6) |
|
412
|
|
|
|
|
413
|
|
|
|
|
414
|
|
|
def test_cape_cin_no_lfc(): |
|
415
|
|
|
"""Tests that CAPE is zero with no LFC.""" |
|
416
|
|
|
p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar |
|
417
|
|
|
temperature = np.array([22.2, 24.6, 22., 20.4, 18., -10.]) * units.celsius |
|
418
|
|
|
dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius |
|
419
|
|
|
parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC') |
|
420
|
|
|
cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof) |
|
421
|
|
|
assert_almost_equal(cape, 0.0 * units('joule / kilogram'), 6) |
|
422
|
|
|
assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 6) |
|
423
|
|
|
|
|
424
|
|
|
|
|
425
|
|
|
def test_find_append_zero_crossings(): |
|
426
|
|
|
"""Tests finding and appending zero crossings of an x, y series.""" |
|
427
|
|
|
x = np.arange(11) * units.hPa |
|
428
|
|
|
y = np.array([3, 2, 1, -1, 2, 2, 0, 1, 0, -1, 2]) * units.degC |
|
429
|
|
|
x2, y2 = _find_append_zero_crossings(x, y) |
|
430
|
|
|
|
|
431
|
|
|
x_truth = np.array([0., 1., 2., 2.5, 3., 3.33333333, 4., 5., |
|
432
|
|
|
6., 7., 8., 9., 9.33333333, 10.]) * units.hPa |
|
433
|
|
|
y_truth = np.array([3, 2, 1, 0, -1, 0, 2, 2, 0, 1, 0, -1, 0, 2]) * units.degC |
|
434
|
|
|
assert_array_almost_equal(x2, x_truth, 6) |
|
435
|
|
|
assert_almost_equal(y2, y_truth, 6) |
|
436
|
|
|
|
|
437
|
|
View Code Duplication |
|
|
|
|
|
|
|
438
|
|
|
def test_most_unstable_parcel(): |
|
439
|
|
|
"""Tests calculating the most unstable parcel.""" |
|
440
|
|
|
levels = np.array([1000., 959., 867.9]) * units.mbar |
|
441
|
|
|
temperatures = np.array([18.2, 22.2, 17.4]) * units.celsius |
|
442
|
|
|
dewpoints = np.array([19., 19., 14.3]) * units.celsius |
|
443
|
|
|
ret = most_unstable_parcel(levels, temperatures, dewpoints, depth=100 * units.hPa) |
|
444
|
|
|
assert_almost_equal(ret[0], 959.0 * units.hPa, 6) |
|
445
|
|
|
assert_almost_equal(ret[1], 22.2 * units.degC, 6) |
|
446
|
|
|
assert_almost_equal(ret[2], 19.0 * units.degC, 6) |
|
447
|
|
|
|
|
448
|
|
|
|
|
449
|
|
|
def test_isentropic_pressure(): |
|
450
|
|
|
"""Test calculation of isentropic pressure function.""" |
|
451
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
452
|
|
|
tmp = np.ones((4, 5, 5)) |
|
453
|
|
View Code Duplication |
tmp[0, :] = 296. |
|
|
|
|
|
|
454
|
|
|
tmp[1, :] = 292. |
|
455
|
|
|
tmp[2, :] = 290 |
|
456
|
|
|
tmp[3, :] = 288. |
|
457
|
|
|
tmpk = tmp * units.kelvin |
|
458
|
|
|
isentlev = [296.] * units.kelvin |
|
459
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk) |
|
460
|
|
|
trueprs = 1000. * units.hPa |
|
461
|
|
|
assert_almost_equal(isentprs[0].shape, (1, 5, 5), 3) |
|
462
|
|
|
assert_almost_equal(isentprs[0], trueprs, 3) |
|
463
|
|
|
|
|
464
|
|
|
|
|
465
|
|
|
def test_isentropic_pressure_p_increase(): |
|
466
|
|
|
"""Test calculation of isentropic pressure function, p increasing order.""" |
|
467
|
|
|
lev = [85000, 90000., 95000., 100000.] * units.Pa |
|
468
|
|
View Code Duplication |
tmp = np.ones((4, 5, 5)) |
|
|
|
|
|
|
469
|
|
|
tmp[0, :] = 288. |
|
470
|
|
|
tmp[1, :] = 290. |
|
471
|
|
|
tmp[2, :] = 292. |
|
472
|
|
|
tmp[3, :] = 296. |
|
473
|
|
|
tmpk = tmp * units.kelvin |
|
474
|
|
|
isentlev = [296.] * units.kelvin |
|
475
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk) |
|
476
|
|
|
trueprs = 1000. * units.hPa |
|
477
|
|
|
assert_almost_equal(isentprs[0], trueprs, 3) |
|
478
|
|
|
|
|
479
|
|
|
|
|
480
|
|
|
def test_isentropic_pressure_adition_args(): |
|
481
|
|
|
"""Test calculation of isentropic pressure function, additional args.""" |
|
482
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
483
|
|
|
tmp = np.ones((4, 5, 5)) |
|
484
|
|
|
tmp[0, :] = 296. |
|
485
|
|
|
tmp[1, :] = 292. |
|
486
|
|
|
tmp[2, :] = 290. |
|
487
|
|
|
tmp[3, :] = 288. |
|
488
|
|
|
rh = np.ones((4, 5, 5)) |
|
489
|
|
View Code Duplication |
rh[0, :] = 100. |
|
|
|
|
|
|
490
|
|
|
rh[1, :] = 80. |
|
491
|
|
|
rh[2, :] = 40. |
|
492
|
|
|
rh[3, :] = 20. |
|
493
|
|
|
relh = rh * units.percent |
|
494
|
|
|
tmpk = tmp * units.kelvin |
|
495
|
|
|
isentlev = [296.] * units.kelvin |
|
496
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh) |
|
497
|
|
|
truerh = 100. * units.percent |
|
498
|
|
|
assert_almost_equal(isentprs[1], truerh, 3) |
|
499
|
|
|
|
|
500
|
|
|
|
|
501
|
|
|
def test_isentropic_pressure_tmp_out(): |
|
502
|
|
|
"""Test calculation of isentropic pressure function, temperature output.""" |
|
503
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
504
|
|
View Code Duplication |
tmp = np.ones((4, 5, 5)) |
|
|
|
|
|
|
505
|
|
|
tmp[0, :] = 296. |
|
506
|
|
|
tmp[1, :] = 292. |
|
507
|
|
|
tmp[2, :] = 290. |
|
508
|
|
|
tmp[3, :] = 288. |
|
509
|
|
|
tmpk = tmp * units.kelvin |
|
510
|
|
|
isentlev = [296.] * units.kelvin |
|
511
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk, tmpk_out=True) |
|
512
|
|
|
truetmp = 296. * units.kelvin |
|
513
|
|
|
assert_almost_equal(isentprs[1], truetmp, 3) |
|
514
|
|
|
|
|
515
|
|
|
|
|
516
|
|
|
def test_isentropic_pressure_p_increase_rh_out(): |
|
517
|
|
|
"""Test calculation of isentropic pressure function, p increasing order.""" |
|
518
|
|
|
lev = [85000., 90000., 95000., 100000.] * units.Pa |
|
519
|
|
|
tmp = np.ones((4, 5, 5)) |
|
520
|
|
|
tmp[0, :] = 288. |
|
521
|
|
|
tmp[1, :] = 290. |
|
522
|
|
|
tmp[2, :] = 292. |
|
523
|
|
|
tmp[3, :] = 296. |
|
524
|
|
|
tmpk = tmp * units.kelvin |
|
525
|
|
View Code Duplication |
rh = np.ones((4, 5, 5)) |
|
|
|
|
|
|
526
|
|
|
rh[0, :] = 20. |
|
527
|
|
|
rh[1, :] = 40. |
|
528
|
|
|
rh[2, :] = 80. |
|
529
|
|
|
rh[3, :] = 100. |
|
530
|
|
|
relh = rh * units.percent |
|
531
|
|
|
isentlev = 296. * units.kelvin |
|
532
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh) |
|
533
|
|
|
truerh = 100. * units.percent |
|
534
|
|
|
assert_almost_equal(isentprs[1], truerh, 3) |
|
535
|
|
|
|
|
536
|
|
|
|
|
537
|
|
|
def test_isentropic_pressure_interp(): |
|
538
|
|
|
"""Test calculation of isentropic pressure function.""" |
|
539
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
540
|
|
View Code Duplication |
tmp = np.ones((4, 5, 5)) |
|
|
|
|
|
|
541
|
|
|
tmp[0, :] = 296. |
|
542
|
|
|
tmp[1, :] = 292. |
|
543
|
|
|
tmp[2, :] = 290 |
|
544
|
|
|
tmp[3, :] = 288. |
|
545
|
|
|
tmpk = tmp * units.kelvin |
|
546
|
|
|
isentlev = [296., 297] * units.kelvin |
|
547
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk) |
|
548
|
|
|
trueprs = 936.18057 * units.hPa |
|
549
|
|
|
assert_almost_equal(isentprs[0][1], trueprs, 3) |
|
550
|
|
|
|
|
551
|
|
|
|
|
552
|
|
|
def test_isentropic_pressure_adition_args_interp(): |
|
553
|
|
|
"""Test calculation of isentropic pressure function, additional args.""" |
|
554
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
555
|
|
|
tmp = np.ones((4, 5, 5)) |
|
556
|
|
|
tmp[0, :] = 296. |
|
557
|
|
|
tmp[1, :] = 292. |
|
558
|
|
|
tmp[2, :] = 290. |
|
559
|
|
|
tmp[3, :] = 288. |
|
560
|
|
|
rh = np.ones((4, 5, 5)) |
|
561
|
|
View Code Duplication |
rh[0, :] = 100. |
|
|
|
|
|
|
562
|
|
|
rh[1, :] = 80. |
|
563
|
|
|
rh[2, :] = 40. |
|
564
|
|
|
rh[3, :] = 20. |
|
565
|
|
|
relh = rh * units.percent |
|
566
|
|
|
tmpk = tmp * units.kelvin |
|
567
|
|
|
isentlev = [296., 297.] * units.kelvin |
|
568
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh) |
|
569
|
|
|
truerh = 69.171 * units.percent |
|
570
|
|
|
assert_almost_equal(isentprs[1][1], truerh, 3) |
|
571
|
|
|
|
|
572
|
|
|
|
|
573
|
|
|
def test_isentropic_pressure_tmp_out_interp(): |
|
574
|
|
|
"""Test calculation of isentropic pressure function, temperature output.""" |
|
575
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
576
|
|
|
tmp = np.ones((4, 5, 5)) |
|
577
|
|
|
tmp[0, :] = 296. |
|
578
|
|
|
tmp[1, :] = 292. |
|
579
|
|
|
tmp[2, :] = 290. |
|
580
|
|
|
tmp[3, :] = 288. |
|
581
|
|
|
tmpk = tmp * units.kelvin |
|
582
|
|
|
isentlev = [296., 297.] * units.kelvin |
|
583
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk, tmpk_out=True) |
|
584
|
|
|
truetmp = 291.4579 * units.kelvin |
|
585
|
|
|
assert_almost_equal(isentprs[1][1], truetmp, 3) |
|
586
|
|
|
|
|
587
|
|
|
|
|
588
|
|
|
def test_isentropic_pressure_data_bounds_error(): |
|
589
|
|
|
"""Test calculation of isentropic pressure function, error for data out of bounds.""" |
|
590
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
591
|
|
|
tmp = np.ones((4, 5, 5)) |
|
592
|
|
|
tmp[0, :] = 296. |
|
593
|
|
|
tmp[1, :] = 292. |
|
594
|
|
|
tmp[2, :] = 290. |
|
595
|
|
|
tmp[3, :] = 288. |
|
596
|
|
|
tmpk = tmp * units.kelvin |
|
597
|
|
|
isentlev = [296., 350.] * units.kelvin |
|
598
|
|
|
with pytest.raises(ValueError): |
|
599
|
|
|
isentropic_interpolation(isentlev, lev, tmpk) |
|
600
|
|
|
|
|
601
|
|
|
|
|
602
|
|
|
def test_isentropic_pressure_4d(): |
|
603
|
|
|
"""Test calculation of isentropic pressure function.""" |
|
604
|
|
|
lev = [100000., 95000., 90000., 85000.] * units.Pa |
|
605
|
|
|
tmp = np.ones((3, 4, 5, 5)) |
|
606
|
|
|
tmp[:, 0, :] = 296. |
|
607
|
|
|
tmp[:, 1, :] = 292. |
|
608
|
|
|
tmp[:, 2, :] = 290 |
|
609
|
|
|
tmp[:, 3, :] = 288. |
|
610
|
|
|
tmpk = tmp * units.kelvin |
|
611
|
|
|
rh = np.ones((3, 4, 5, 5)) |
|
612
|
|
|
rh[:, 0, :] = 100. |
|
613
|
|
|
rh[:, 1, :] = 80. |
|
614
|
|
|
rh[:, 2, :] = 40. |
|
615
|
|
|
rh[:, 3, :] = 20. |
|
616
|
|
|
relh = rh * units.percent |
|
617
|
|
|
isentlev = [296., 297., 300.] * units.kelvin |
|
618
|
|
|
isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh, axis=1) |
|
619
|
|
|
trueprs = 1000. * units.hPa |
|
620
|
|
|
trueprs2 = 936.18057 * units.hPa |
|
621
|
|
|
trueprs3 = 879.446 * units.hPa |
|
622
|
|
|
truerh = 69.171 * units.percent |
|
623
|
|
|
assert_almost_equal(isentprs[0].shape, (3, 3, 5, 5), 3) |
|
624
|
|
|
assert_almost_equal(isentprs[0][:, 0, :], trueprs, 3) |
|
625
|
|
|
assert_almost_equal(isentprs[0][:, 1, :], trueprs2, 3) |
|
626
|
|
|
assert_almost_equal(isentprs[0][:, 2, :], trueprs3, 3) |
|
627
|
|
|
assert_almost_equal(isentprs[1][:, 1, ], truerh, 3) |
|
628
|
|
|
|
Unidata /
MetPy
Loading history...
