|
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 (density, dewpoint, dewpoint_rh, dry_lapse, el, |
|
10
|
|
|
equivalent_potential_temperature, lcl, lfc, mixing_ratio, moist_lapse, |
|
11
|
|
|
parcel_profile, potential_temperature, |
|
12
|
|
|
psychrometric_vapor_pressure_wet, relative_humidity_wet_psychrometric, |
|
13
|
|
|
saturation_mixing_ratio, saturation_vapor_pressure, vapor_pressure, |
|
14
|
|
|
virtual_potential_temperature, virtual_temperature) |
|
15
|
|
|
|
|
16
|
|
|
from metpy.testing import assert_almost_equal, assert_array_almost_equal |
|
17
|
|
|
from metpy.units import units |
|
18
|
|
|
|
|
19
|
|
|
|
|
20
|
|
|
def assertIsUnitNan(value, units): |
|
21
|
|
|
"""Helper to check for nan with proper units.""" |
|
22
|
|
|
if not np.isnan(value): |
|
23
|
|
|
pytest.fail("{} is not np.nan".format(value)) |
|
24
|
|
|
if not hasattr(value, 'units'): |
|
25
|
|
|
pytest.fail("No units attached") |
|
26
|
|
|
if not value.units == units: |
|
27
|
|
|
pytest.fail("{} units are not {}".format(value.units, units)) |
|
28
|
|
|
return True |
|
29
|
|
|
|
|
30
|
|
|
def test_potential_temperature(): |
|
31
|
|
|
"""Test potential_temperature calculation.""" |
|
32
|
|
|
temp = np.array([278., 283., 291., 298.]) * units.kelvin |
|
33
|
|
|
pres = np.array([900., 500., 300., 100.]) * units.mbar |
|
34
|
|
|
real_th = np.array([286.493, 344.961, 410.4335, 575.236]) * units.kelvin |
|
35
|
|
|
assert_array_almost_equal(potential_temperature(pres, temp), real_th, 3) |
|
36
|
|
|
|
|
37
|
|
|
|
|
38
|
|
|
def test_scalar(): |
|
39
|
|
|
"""Test potential_temperature accepts scalar values.""" |
|
40
|
|
|
assert_almost_equal(potential_temperature(1000. * units.mbar, 293. * units.kelvin), |
|
41
|
|
|
293. * units.kelvin, 4) |
|
42
|
|
|
assert_almost_equal(potential_temperature(800. * units.mbar, 293. * units.kelvin), |
|
43
|
|
|
312.2828 * units.kelvin, 4) |
|
44
|
|
|
|
|
45
|
|
|
|
|
46
|
|
|
def test_fahrenheit(): |
|
47
|
|
|
"""Test that potential_temperature handles temperature values in Fahrenheit.""" |
|
48
|
|
|
assert_almost_equal(potential_temperature(800. * units.mbar, 68. * units.degF), |
|
49
|
|
|
(312.444 * units.kelvin).to(units.degF), 2) |
|
50
|
|
|
|
|
51
|
|
|
|
|
52
|
|
|
def test_pot_temp_inhg(): |
|
53
|
|
|
"""Test that potential_temperature can handle pressure not in mb (issue #165).""" |
|
54
|
|
|
assert_almost_equal(potential_temperature(29.92 * units.inHg, 29 * units.degC), |
|
55
|
|
|
301.019735 * units.kelvin, 4) |
|
56
|
|
|
|
|
57
|
|
|
|
|
58
|
|
|
def test_dry_lapse(): |
|
59
|
|
|
"""Test dry_lapse calculation.""" |
|
60
|
|
|
levels = np.array([1000, 900, 864.89]) * units.mbar |
|
61
|
|
|
temps = dry_lapse(levels, 303.15 * units.kelvin) |
|
62
|
|
|
assert_array_almost_equal(temps, |
|
63
|
|
|
np.array([303.15, 294.16, 290.83]) * units.kelvin, 2) |
|
64
|
|
|
|
|
65
|
|
|
|
|
66
|
|
|
def test_dry_lapse_2_levels(): |
|
67
|
|
|
"""Test dry_lapse calculation when given only two levels.""" |
|
68
|
|
|
temps = dry_lapse(np.array([1000., 500.]) * units.mbar, 293. * units.kelvin) |
|
69
|
|
|
assert_array_almost_equal(temps, [293., 240.3723] * units.kelvin, 4) |
|
70
|
|
|
|
|
71
|
|
|
|
|
72
|
|
|
def test_moist_lapse(): |
|
73
|
|
|
"""Test moist_lapse calculation.""" |
|
74
|
|
|
temp = moist_lapse(np.array([1000., 800., 600., 500., 400.]) * units.mbar, |
|
75
|
|
|
293. * units.kelvin) |
|
76
|
|
|
true_temp = np.array([293, 284.64, 272.81, 264.42, 252.91]) * units.kelvin |
|
77
|
|
|
assert_array_almost_equal(temp, true_temp, 2) |
|
78
|
|
|
|
|
79
|
|
|
|
|
80
|
|
|
def test_moist_lapse_degc(): |
|
81
|
|
|
"""Test moist_lapse with Celsius temperatures.""" |
|
82
|
|
|
temp = moist_lapse(np.array([1000., 800., 600., 500., 400.]) * units.mbar, |
|
83
|
|
|
19.85 * units.degC) |
|
84
|
|
|
true_temp = np.array([293, 284.64, 272.81, 264.42, 252.91]) * units.kelvin |
|
85
|
|
|
assert_array_almost_equal(temp, true_temp, 2) |
|
86
|
|
|
|
|
87
|
|
|
|
|
88
|
|
|
def test_parcel_profile(): |
|
89
|
|
|
"""Test parcel profile calculation.""" |
|
90
|
|
|
levels = np.array([1000., 900., 800., 700., 600., 500., 400.]) * units.mbar |
|
91
|
|
|
true_prof = np.array([303.15, 294.16, 288.026, 283.073, 277.058, 269.402, |
|
92
|
|
|
258.966]) * units.kelvin |
|
93
|
|
|
|
|
94
|
|
|
prof = parcel_profile(levels, 30. * units.degC, 20. * units.degC) |
|
95
|
|
|
assert_array_almost_equal(prof, true_prof, 2) |
|
96
|
|
|
|
|
97
|
|
|
|
|
98
|
|
|
def test_parcel_profile_saturated(): |
|
99
|
|
|
"""Test parcel_profile works when LCL in levels (issue #232).""" |
|
100
|
|
|
levels = np.array([1000., 700., 500.]) * units.mbar |
|
101
|
|
|
true_prof = np.array([296.95, 284.381, 271.123]) * units.kelvin |
|
102
|
|
|
|
|
103
|
|
|
prof = parcel_profile(levels, 23.8 * units.degC, 23.8 * units.degC) |
|
104
|
|
|
assert_array_almost_equal(prof, true_prof, 2) |
|
105
|
|
|
|
|
106
|
|
|
|
|
107
|
|
|
def test_sat_vapor_pressure(): |
|
108
|
|
|
"""Test saturation_vapor_pressure calculation.""" |
|
109
|
|
|
temp = np.array([5., 10., 18., 25.]) * units.degC |
|
110
|
|
|
real_es = np.array([8.72, 12.27, 20.63, 31.67]) * units.mbar |
|
111
|
|
|
assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 2) |
|
112
|
|
|
|
|
113
|
|
|
|
|
114
|
|
|
def test_sat_vapor_pressure_scalar(): |
|
115
|
|
|
"""Test saturation_vapor_pressure handles scalar values.""" |
|
116
|
|
|
es = saturation_vapor_pressure(0 * units.degC) |
|
117
|
|
|
assert_almost_equal(es, 6.112 * units.mbar, 3) |
|
118
|
|
|
|
|
119
|
|
|
|
|
120
|
|
|
def test_sat_vapor_pressure_fahrenheit(): |
|
121
|
|
|
"""Test saturation_vapor_pressure handles temperature in Fahrenheit.""" |
|
122
|
|
|
temp = np.array([50., 68.]) * units.degF |
|
123
|
|
|
real_es = np.array([12.2717, 23.3695]) * units.mbar |
|
124
|
|
|
assert_array_almost_equal(saturation_vapor_pressure(temp), real_es, 4) |
|
125
|
|
|
|
|
126
|
|
|
|
|
127
|
|
|
def test_basic_dewpoint_rh(): |
|
128
|
|
|
"""Test dewpoint_rh function.""" |
|
129
|
|
|
temp = np.array([30., 25., 10., 20., 25.]) * units.degC |
|
130
|
|
|
rh = np.array([30., 45., 55., 80., 85.]) / 100. |
|
131
|
|
|
|
|
132
|
|
|
real_td = np.array([11, 12, 1, 16, 22]) * units.degC |
|
133
|
|
|
assert_array_almost_equal(real_td, dewpoint_rh(temp, rh), 0) |
|
134
|
|
|
|
|
135
|
|
|
|
|
136
|
|
|
def test_scalar_dewpoint_rh(): |
|
137
|
|
|
"""Test dewpoint_rh with scalar values.""" |
|
138
|
|
|
td = dewpoint_rh(10.6 * units.degC, 0.37) |
|
139
|
|
|
assert_almost_equal(td, 26. * units.degF, 0) |
|
140
|
|
|
|
|
141
|
|
|
|
|
142
|
|
|
def test_dewpoint(): |
|
143
|
|
|
"""Test dewpoint calculation.""" |
|
144
|
|
|
assert_almost_equal(dewpoint(6.112 * units.mbar), 0. * units.degC, 2) |
|
145
|
|
|
|
|
146
|
|
|
|
|
147
|
|
|
def test_dewpoint_weird_units(): |
|
148
|
|
|
"""Test dewpoint using non-standard units. |
|
149
|
|
|
|
|
150
|
|
|
Revealed from odd dimensionless units and ending up using numpy.ma math |
|
151
|
|
|
functions instead of numpy ones. |
|
152
|
|
|
""" |
|
153
|
|
|
assert_almost_equal(dewpoint(15825.6 * units('g * mbar / kg')), |
|
154
|
|
|
13.8564 * units.degC, 4) |
|
155
|
|
|
|
|
156
|
|
|
|
|
157
|
|
|
def test_mixing_ratio(): |
|
158
|
|
|
"""Test mixing ratio calculation.""" |
|
159
|
|
|
p = 998. * units.mbar |
|
160
|
|
|
e = 73.75 * units.mbar |
|
161
|
|
|
assert_almost_equal(mixing_ratio(e, p), 0.04963, 2) |
|
162
|
|
|
|
|
163
|
|
|
|
|
164
|
|
|
def test_vapor_pressure(): |
|
165
|
|
|
"""Test vapor pressure calculation.""" |
|
166
|
|
|
assert_almost_equal(vapor_pressure(998. * units.mbar, 0.04963), |
|
167
|
|
|
73.74925 * units.mbar, 5) |
|
168
|
|
|
|
|
169
|
|
|
|
|
170
|
|
|
def test_lcl(): |
|
171
|
|
|
"""Test LCL calculation.""" |
|
172
|
|
|
lcl_pressure, lcl_temperature = lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC) |
|
173
|
|
|
assert_almost_equal(lcl_pressure, 864.761 * units.mbar, 2) |
|
174
|
|
|
assert_almost_equal(lcl_temperature, 17.676 * units.degC, 2) |
|
175
|
|
|
|
|
176
|
|
|
|
|
177
|
|
|
def test_lcl_convergence(): |
|
178
|
|
|
"""Test LCL calculation convergence failure.""" |
|
179
|
|
|
with pytest.raises(RuntimeError): |
|
180
|
|
|
lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC, max_iters=2) |
|
181
|
|
|
|
|
182
|
|
|
|
|
183
|
|
|
def test_lfc_basic(): |
|
184
|
|
|
"""Test LFC calculation.""" |
|
185
|
|
|
levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar |
|
186
|
|
|
temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -49.]) * units.celsius |
|
187
|
|
|
dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius |
|
188
|
|
|
l = lfc(levels, temperatures, dewpoints) |
|
189
|
|
|
assert_almost_equal(l[0], 727.468 * units.mbar, 2) |
|
190
|
|
|
assert_almost_equal(l[1], 9.705 * units.celsius, 2) |
|
191
|
|
|
|
|
192
|
|
|
|
|
193
|
|
|
def test_no_lfc(): |
|
194
|
|
|
"""Test LFC calculation when there is no LFC in the data.""" |
|
195
|
|
|
levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar |
|
196
|
|
|
temperatures = np.array([22.2, 17.4, 14.6, 1.4, -17.6, -39.4, -52.5]) * units.celsius |
|
197
|
|
|
dewpoints = np.array([9., 4.3, -21.2, -26.7, -31., -53.3, -66.7]) * units.celsius |
|
198
|
|
|
lfc_pressure, lfc_temperature = lfc(levels, temperatures, dewpoints) |
|
199
|
|
|
assert assertIsUnitNan(lfc_pressure, levels.units) |
|
200
|
|
|
assert assertIsUnitNan(lfc_temperature, temperatures.units) |
|
201
|
|
|
|
|
202
|
|
|
|
|
203
|
|
|
|
|
204
|
|
|
def test_lfc_inversion(): |
|
205
|
|
|
"""Test LFC when there is an inversion to be sure we don't pick that.""" |
|
206
|
|
|
levels = np.array([963., 789., 782.3, 754.8, 728.1, 727., 700., |
|
207
|
|
|
571., 450., 300., 248.]) * units.mbar |
|
208
|
|
|
temperatures = np.array([25.4, 18.4, 17.8, 15.4, 12.9, 12.8, |
|
209
|
|
|
10., -3.9, -16.3, -41.1, -51.5]) * units.celsius |
|
210
|
|
|
dewpoints = np.array([20.4, 0.4, -0.5, -4.3, -8., -8.2, -9., |
|
211
|
|
|
-23.9, -33.3, -54.1, -63.5]) * units.celsius |
|
212
|
|
|
l = lfc(levels, temperatures, dewpoints) |
|
213
|
|
|
assert_almost_equal(l[0], 706.0103 * units.mbar, 2) |
|
214
|
|
|
assert_almost_equal(l[1], 10.6232 * units.celsius, 2) |
|
215
|
|
|
|
|
216
|
|
|
|
|
217
|
|
|
def test_saturation_mixing_ratio(): |
|
218
|
|
|
"""Test saturation mixing ratio calculation.""" |
|
219
|
|
|
p = 999. * units.mbar |
|
220
|
|
|
t = 288. * units.kelvin |
|
221
|
|
|
assert_almost_equal(saturation_mixing_ratio(p, t), .01068, 3) |
|
222
|
|
|
|
|
223
|
|
|
|
|
224
|
|
|
def test_equivalent_potential_temperature(): |
|
225
|
|
|
"""Test equivalent potential temperature calculation.""" |
|
226
|
|
|
p = 999. * units.mbar |
|
227
|
|
|
t = 288. * units.kelvin |
|
228
|
|
|
ept = equivalent_potential_temperature(p, t) |
|
229
|
|
|
assert_almost_equal(ept, 315.9548 * units.kelvin, 3) |
|
230
|
|
|
|
|
231
|
|
|
|
|
232
|
|
|
def test_virtual_temperature(): |
|
233
|
|
|
"""Test virtual temperature calculation.""" |
|
234
|
|
|
t = 288. * units.kelvin |
|
235
|
|
|
qv = .0016 # kg/kg |
|
236
|
|
|
tv = virtual_temperature(t, qv) |
|
237
|
|
|
assert_almost_equal(tv, 288.2796 * units.kelvin, 3) |
|
238
|
|
|
|
|
239
|
|
|
|
|
240
|
|
|
def test_virtual_potential_temperature(): |
|
241
|
|
|
"""Test virtual potential temperature calculation.""" |
|
242
|
|
|
p = 999. * units.mbar |
|
243
|
|
|
t = 288. * units.kelvin |
|
244
|
|
|
qv = .0016 # kg/kg |
|
245
|
|
|
theta_v = virtual_potential_temperature(p, t, qv) |
|
246
|
|
|
assert_almost_equal(theta_v, 288.3620 * units.kelvin, 3) |
|
247
|
|
|
|
|
248
|
|
|
|
|
249
|
|
|
def test_density(): |
|
250
|
|
|
"""Test density calculation.""" |
|
251
|
|
|
p = 999. * units.mbar |
|
252
|
|
|
t = 288. * units.kelvin |
|
253
|
|
|
qv = .0016 # kg/kg |
|
254
|
|
|
rho = density(p, t, qv).to(units.kilogram / units.meter ** 3) |
|
255
|
|
|
assert_almost_equal(rho, 1.2072 * (units.kilogram / units.meter ** 3), 3) |
|
256
|
|
|
|
|
257
|
|
|
|
|
258
|
|
|
def test_el(): |
|
259
|
|
|
"""Test equilibrium layer calculation.""" |
|
260
|
|
|
levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar |
|
261
|
|
|
temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius |
|
262
|
|
|
dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius |
|
263
|
|
|
el_pressure, el_temperature = el(levels, temperatures, dewpoints) |
|
264
|
|
|
assert_almost_equal(el_pressure, 520.8420 * units.mbar, 3) |
|
265
|
|
|
assert_almost_equal(el_temperature, -11.7055 * units.degC, 3) |
|
266
|
|
|
|
|
267
|
|
|
|
|
268
|
|
|
def test_no_el(): |
|
269
|
|
|
"""Test equilibrium layer calculation when there is no EL in the data.""" |
|
270
|
|
|
levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar |
|
271
|
|
|
temperatures = np.array([22.2, 17.4, 14.6, 1.4, -17.6, -39.4, -52.5]) * units.celsius |
|
272
|
|
|
dewpoints = np.array([19., 14.3, -11.2, -16.7, -21., -43.3, -56.7]) * units.celsius |
|
273
|
|
|
el_pressure, el_temperature = el(levels, temperatures, dewpoints) |
|
274
|
|
|
assert assertIsUnitNan(el_pressure, levels.units) |
|
275
|
|
|
assert assertIsUnitNan(el_temperature, temperatures.units) |
|
276
|
|
|
|
|
277
|
|
|
|
|
278
|
|
|
def test_wet_psychrometric_vapor_pressure(): |
|
279
|
|
|
"""Test calculation of vapor pressure from wet and dry bulb temperatures.""" |
|
280
|
|
|
p = 1013.25 * units.mbar |
|
281
|
|
|
dry_bulb_temperature = 20. * units.degC |
|
282
|
|
|
wet_bulb_temperature = 18. * units.degC |
|
283
|
|
|
psychrometric_vapor_pressure = psychrometric_vapor_pressure_wet(dry_bulb_temperature, |
|
284
|
|
|
wet_bulb_temperature, p) |
|
285
|
|
|
assert_almost_equal(psychrometric_vapor_pressure, 19.3673 * units.mbar, 3) |
|
286
|
|
|
|
|
287
|
|
|
|
|
288
|
|
|
def test_wet_psychrometric_rh(): |
|
289
|
|
|
"""Test calculation of relative humidity from wet and dry bulb temperatures.""" |
|
290
|
|
|
p = 1013.25 * units.mbar |
|
291
|
|
|
dry_bulb_temperature = 20. * units.degC |
|
292
|
|
|
wet_bulb_temperature = 18. * units.degC |
|
293
|
|
|
psychrometric_rh = relative_humidity_wet_psychrometric(dry_bulb_temperature, |
|
294
|
|
|
wet_bulb_temperature, p) |
|
295
|
|
|
assert_almost_equal(psychrometric_rh, 82.8747 * units.percent, 3) |
|
296
|
|
|
|
|
297
|
|
|
|
|
298
|
|
|
def test_wet_psychrometric_rh_kwargs(): |
|
299
|
|
|
"""Test calculation of relative humidity from wet and dry bulb temperatures.""" |
|
300
|
|
|
p = 1013.25 * units.mbar |
|
301
|
|
|
dry_bulb_temperature = 20. * units.degC |
|
302
|
|
|
wet_bulb_temperature = 18. * units.degC |
|
303
|
|
|
coeff = 6.1e-4 / units.kelvin |
|
304
|
|
|
psychrometric_rh = relative_humidity_wet_psychrometric(dry_bulb_temperature, |
|
305
|
|
|
wet_bulb_temperature, p, |
|
306
|
|
|
psychrometer_coefficient=coeff) |
|
307
|
|
|
assert_almost_equal(psychrometric_rh, 82.9701 * units.percent, 3) |
|
308
|
|
|
|
Unidata /
MetPy
