Completed
Push — master ( e1d789...31695a )
by Ryan
18s
created

test_el_small_surface_instability()   A

Complexity

Conditions 3

Size

Total Lines 14

Duplication

Lines 0
Ratio 0 %

Importance

Changes 1
Bugs 0 Features 0
Metric Value
cc 3
dl 0
loc 14
rs 9.4285
c 1
b 0
f 0
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, mixed_layer, mixed_parcel, mixing_ratio,
13
                        mixing_ratio_from_specific_humidity, moist_lapse,
14
                        most_unstable_cape_cin, most_unstable_parcel,
15
                        parcel_profile, potential_temperature,
16
                        psychrometric_vapor_pressure_wet,
17
                        relative_humidity_from_mixing_ratio,
18
                        relative_humidity_from_specific_humidity,
19
                        relative_humidity_wet_psychrometric,
20
                        saturation_mixing_ratio,
21
                        saturation_vapor_pressure,
22
                        surface_based_cape_cin, vapor_pressure,
23
                        virtual_potential_temperature, virtual_temperature)
24
25
from metpy.calc.thermo import _find_append_zero_crossings
26
from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan
27
from metpy.units import units
28
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_percent_dewpoint_rh():
143
    """Test dewpoint_rh with rh in percent."""
144
    td = dewpoint_rh(10.6 * units.degC, 37 * units.percent)
145
    assert_almost_equal(td, 26. * units.degF, 0)
146
147
148
def test_warning_dewpoint_rh():
149
    """Test that warning is raised for >120% RH."""
150
    with pytest.warns(UserWarning):
151
        dewpoint_rh(10.6 * units.degC, 50)
152
153
154
def test_dewpoint():
155
    """Test dewpoint calculation."""
156
    assert_almost_equal(dewpoint(6.112 * units.mbar), 0. * units.degC, 2)
157
158
159
def test_dewpoint_weird_units():
160
    """Test dewpoint using non-standard units.
161
162
    Revealed from odd dimensionless units and ending up using numpy.ma math
163
    functions instead of numpy ones.
164
    """
165
    assert_almost_equal(dewpoint(15825.6 * units('g * mbar / kg')),
166
                        13.8564 * units.degC, 4)
167
168
169
def test_mixing_ratio():
170
    """Test mixing ratio calculation."""
171
    p = 998. * units.mbar
172
    e = 73.75 * units.mbar
173
    assert_almost_equal(mixing_ratio(e, p), 0.04963, 2)
174
175
176
def test_vapor_pressure():
177
    """Test vapor pressure calculation."""
178
    assert_almost_equal(vapor_pressure(998. * units.mbar, 0.04963),
179 View Code Duplication
                        73.74925 * units.mbar, 5)
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
180
181
182
def test_lcl():
183
    """Test LCL calculation."""
184
    lcl_pressure, lcl_temperature = lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC)
185
    assert_almost_equal(lcl_pressure, 864.761 * units.mbar, 2)
186
    assert_almost_equal(lcl_temperature, 17.676 * units.degC, 2)
187
188
189
def test_lcl_convergence():
190
    """Test LCL calculation convergence failure."""
191
    with pytest.raises(RuntimeError):
192
        lcl(1000. * units.mbar, 30. * units.degC, 20. * units.degC, max_iters=2)
193
194
195
def test_lfc_basic():
196
    """Test LFC calculation."""
197
    levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
198
    temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -49.]) * units.celsius
199
    dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
200
    l = lfc(levels, temperatures, dewpoints)
201
    assert_almost_equal(l[0], 727.468 * units.mbar, 2)
202
    assert_almost_equal(l[1], 9.705 * units.celsius, 2)
203
204
205
def test_no_lfc():
206
    """Test LFC calculation when there is no LFC in the data."""
207
    levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar
208
    temperatures = np.array([22.2, 17.4, 14.6, 1.4, -17.6, -39.4, -52.5]) * units.celsius
209
    dewpoints = np.array([9., 4.3, -21.2, -26.7, -31., -53.3, -66.7]) * units.celsius
210
    lfc_pressure, lfc_temperature = lfc(levels, temperatures, dewpoints)
211
    assert assert_nan(lfc_pressure, levels.units)
212
    assert assert_nan(lfc_temperature, temperatures.units)
213
214
215
def test_lfc_inversion():
216
    """Test LFC when there is an inversion to be sure we don't pick that."""
217
    levels = np.array([963., 789., 782.3, 754.8, 728.1, 727., 700.,
218
                       571., 450., 300., 248.]) * units.mbar
219
    temperatures = np.array([25.4, 18.4, 17.8, 15.4, 12.9, 12.8,
220
                             10., -3.9, -16.3, -41.1, -51.5]) * units.celsius
221
    dewpoints = np.array([20.4, 0.4, -0.5, -4.3, -8., -8.2, -9.,
222
                          -23.9, -33.3, -54.1, -63.5]) * units.celsius
223
    l = lfc(levels, temperatures, dewpoints)
224
    assert_almost_equal(l[0], 706.0103 * units.mbar, 2)
225
    assert_almost_equal(l[1], 10.6232 * units.celsius, 2)
226
227
228
def test_lfc_equals_lcl():
229
    """Test LFC when there is no cap and the lfc is equal to the lcl."""
230
    levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1,
231
                       748., 732.2, 700., 654.8]) * units.mbar
232
    temperatures = np.array([29.4, 28.7, 25.2, 22.4, 19.4, 16.8,
233
                             14.3, 13.2, 12.6, 11.4, 7.1]) * units.celsius
234
    dewpoints = np.array([18.4, 18.1, 16.6, 15.4, 13.2, 11.4, 9.6,
235
                          8.8, 0., -18.6, -22.9]) * units.celsius
236
    l = lfc(levels, temperatures, dewpoints)
237
    assert_almost_equal(l[0], 777.0333 * units.mbar, 2)
238
    assert_almost_equal(l[1], 15.8714 * units.celsius, 2)
239
240
241
def test_lfc_sfc_precision():
242
    """Test LFC when there are precision issues with the parcel path."""
243
    levels = np.array([839., 819.4, 816., 807., 790.7, 763., 736.2,
244
                       722., 710.1, 700.]) * units.mbar
245
    temperatures = np.array([20.6, 22.3, 22.6, 22.2, 20.9, 18.7, 16.4,
246
                             15.2, 13.9, 12.8]) * units.celsius
247
    dewpoints = np.array([10.6, 8., 7.6, 6.2, 5.7, 4.7, 3.7, 3.2, 3., 2.8]) * units.celsius
248
    l = lfc(levels, temperatures, dewpoints)
249
    assert assert_nan(l[0], levels.units)
250
    assert assert_nan(l[1], temperatures.units)
251
252
253
def test_saturation_mixing_ratio():
254
    """Test saturation mixing ratio calculation."""
255
    p = 999. * units.mbar
256
    t = 288. * units.kelvin
257
    assert_almost_equal(saturation_mixing_ratio(p, t), .01068, 3)
258
259
260
def test_equivalent_potential_temperature():
261
    """Test equivalent potential temperature calculation."""
262
    p = 1000 * units.mbar
263
    t = 293. * units.kelvin
264
    td = 280. * units.kelvin
265
    ept = equivalent_potential_temperature(p, t, td)
266
    assert_almost_equal(ept, 311.18586467284007 * units.kelvin, 3)
267
268
269
def test_virtual_temperature():
270
    """Test virtual temperature calculation."""
271
    t = 288. * units.kelvin
272
    qv = .0016  # kg/kg
273
    tv = virtual_temperature(t, qv)
274
    assert_almost_equal(tv, 288.2796 * units.kelvin, 3)
275
276
277
def test_virtual_potential_temperature():
278
    """Test virtual potential temperature calculation."""
279
    p = 999. * units.mbar
280
    t = 288. * units.kelvin
281
    qv = .0016  # kg/kg
282
    theta_v = virtual_potential_temperature(p, t, qv)
283
    assert_almost_equal(theta_v, 288.3620 * units.kelvin, 3)
284
285
286
def test_density():
287
    """Test density calculation."""
288
    p = 999. * units.mbar
289
    t = 288. * units.kelvin
290
    qv = .0016  # kg/kg
291
    rho = density(p, t, qv).to(units.kilogram / units.meter ** 3)
292
    assert_almost_equal(rho, 1.2072 * (units.kilogram / units.meter ** 3), 3)
293
294
295
def test_el():
296
    """Test equilibrium layer calculation."""
297
    levels = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
298
    temperatures = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
299
    dewpoints = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
300
    el_pressure, el_temperature = el(levels, temperatures, dewpoints)
301
    assert_almost_equal(el_pressure, 520.8700 * units.mbar, 3)
302
    assert_almost_equal(el_temperature, -11.7027 * units.degC, 3)
303
304
305
def test_no_el():
306
    """Test equilibrium layer calculation when there is no EL in the data."""
307
    levels = np.array([959., 867.9, 779.2, 647.5, 472.5, 321.9, 251.]) * units.mbar
308
    temperatures = np.array([22.2, 17.4, 14.6, 1.4, -17.6, -39.4, -52.5]) * units.celsius
309
    dewpoints = np.array([19., 14.3, -11.2, -16.7, -21., -43.3, -56.7]) * units.celsius
310
    el_pressure, el_temperature = el(levels, temperatures, dewpoints)
311
    assert assert_nan(el_pressure, levels.units)
312
    assert assert_nan(el_temperature, temperatures.units)
313
314
315
def test_no_el_multi_crossing():
316
    """Test el calculation with no el and severel parcel path-profile crossings."""
317
    levels = np.array([918., 911., 880., 873.9, 850., 848., 843.5, 818., 813.8, 785.,
318
                       773., 763., 757.5, 730.5, 700., 679., 654.4, 645.,
319
                       643.9]) * units.mbar
320
    temperatures = np.array([24.2, 22.8, 19.6, 19.1, 17., 16.8, 16.5, 15., 14.9, 14.4, 16.4,
321
                             16.2, 15.7, 13.4, 10.6, 8.4, 5.7, 4.6, 4.5]) * units.celsius
322
    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,
323
                          3.2, 2.6, -0.6, -4.4, -6.6, -9.3, -10.4, -10.5]) * units.celsius
324
    el_pressure, el_temperature = el(levels, temperatures, dewpoints)
325
    assert assert_nan(el_pressure, levels.units)
326
    assert assert_nan(el_temperature, temperatures.units)
327
328
329
def test_el_lfc_equals_lcl():
330
    """Test equilibrium layer calculation when the lfc equals the lcl."""
331
    levels = np.array([912., 905.3, 874.4, 850., 815.1, 786.6, 759.1, 748.,
332
                       732.3, 700., 654.8, 606.8, 562.4, 501.8, 500., 482.,
333
                       400., 393.3, 317.1, 307., 300., 252.7, 250., 200.,
334
                       199.3, 197., 190., 172., 156.6, 150., 122.9, 112.,
335
                       106.2, 100.]) * units.mbar
336
    temperatures = np.array([29.4, 28.7, 25.2, 22.4, 19.4, 16.8, 14.3,
337
                             13.2, 12.6, 11.4, 7.1, 2.2, -2.7, -10.1,
338
                             -10.3, -12.4, -23.3, -24.4, -38., -40.1, -41.1,
339
                             -49.8, -50.3, -59.1, -59.1, -59.3, -59.7, -56.3,
340
                             -56.9, -57.1, -59.1, -60.1, -58.6, -56.9]) * units.celsius
341
    dewpoints = np.array([18.4, 18.1, 16.6, 15.4, 13.2, 11.4, 9.6, 8.8, 0.,
342
                          -18.6, -22.9, -27.8, -32.7, -40.1, -40.3, -42.4, -53.3,
343
                          -54.4, -68., -70.1, -70., -70., -70., -70., -70., -70.,
344
                          -70., -70., -70., -70., -70., -70., -70., -70.]) * units.celsius
345
    el_pressure, el_temperature = el(levels, temperatures, dewpoints)
346
    assert_almost_equal(el_pressure, 175.8684 * units.mbar, 3)
347
    assert_almost_equal(el_temperature, -57.0307 * units.degC, 3)
348
349
350
def test_el_small_surface_instability():
351
    """Test that no EL is found when there is a small pocket of instability at the sfc."""
352
    levels = np.array([959., 931.3, 925., 899.3, 892., 867.9, 850., 814.,
353
                       807.9, 790., 779.2, 751.3, 724.3, 700., 655., 647.5,
354
                       599.4, 554.7, 550., 500.]) * units.mbar
355
    temperatures = np.array([22.2, 20.2, 19.8, 18.4, 18., 17.4, 17., 15.4, 15.4,
356
                             15.6, 14.6, 12., 9.4, 7., 2.2, 1.4, -4.2, -9.7,
357
                             -10.3, -14.9]) * units.degC
358
    dewpoints = np.array([20., 18.5, 18.1, 17.9, 17.8, 15.3, 13.5, 6.4, 2.2,
359
                          -10.4, -10.2, -9.8, -9.4, -9., -15.8, -15.7, -14.8, -14.,
360
                          -13.9, -17.9]) * units.degC
361
    el_pressure, el_temperature = el(levels, temperatures, dewpoints)
362
    assert assert_nan(el_pressure, levels.units)
363
    assert assert_nan(el_temperature, temperatures.units)
364 View Code Duplication
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
365
366
def test_no_el_parcel_colder():
367
    """Tests no EL when parcel stays colder than environment. INL 20170925-12Z."""
368
    levels = np.array([974., 946., 925., 877.2, 866., 850., 814.6, 785.,
369
                       756.6, 739., 729.1, 700., 686., 671., 641., 613.,
370
                       603., 586., 571., 559.3, 539., 533., 500., 491.,
371
                       477.9, 413., 390., 378., 345., 336.]) * units.mbar
372
    temperatures = np.array([10., 8.4, 7.6, 5.9, 7.2, 7.6, 6.8, 7.1, 7.7,
373
                             7.8, 7.7, 5.6, 4.6, 3.4, 0.6, -0.9, -1.1, -3.1,
374
                             -4.7, -4.7, -6.9, -7.5, -11.1, -10.9, -12.1, -20.5, -23.5,
375
                             -24.7, -30.5, -31.7]) * units.celsius
376
    dewpoints = np.array([8.9, 8.4, 7.6, 5.9, 7.2, 7., 5., 3.6, 0.3,
377
                          -4.2, -12.8, -12.4, -8.4, -8.6, -6.4, -7.9, -11.1, -14.1,
378
                          -8.8, -28.1, -18.9, -14.5, -15.2, -15.1, -21.6, -41.5, -45.5,
379
                          -29.6, -30.6, -32.1]) * units.celsius
380
    el_pressure, el_temperature = el(levels, temperatures, dewpoints)
381
    assert assert_nan(el_pressure, levels.units)
382
    assert assert_nan(el_temperature, temperatures.units)
383
384
385
def test_wet_psychrometric_vapor_pressure():
386 View Code Duplication
    """Test calculation of vapor pressure from wet and dry bulb temperatures."""
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
387
    p = 1013.25 * units.mbar
388
    dry_bulb_temperature = 20. * units.degC
389
    wet_bulb_temperature = 18. * units.degC
390
    psychrometric_vapor_pressure = psychrometric_vapor_pressure_wet(dry_bulb_temperature,
391
                                                                    wet_bulb_temperature, p)
392
    assert_almost_equal(psychrometric_vapor_pressure, 19.3673 * units.mbar, 3)
393
394
395
def test_wet_psychrometric_rh():
396
    """Test calculation of relative humidity from wet and dry bulb temperatures."""
397
    p = 1013.25 * units.mbar
398
    dry_bulb_temperature = 20. * units.degC
399
    wet_bulb_temperature = 18. * units.degC
400
    psychrometric_rh = relative_humidity_wet_psychrometric(dry_bulb_temperature,
401
                                                           wet_bulb_temperature, p)
402
    assert_almost_equal(psychrometric_rh, 82.8747 * units.percent, 3)
403
404
405
def test_wet_psychrometric_rh_kwargs():
406
    """Test calculation of relative humidity from wet and dry bulb temperatures."""
407
    p = 1013.25 * units.mbar
408
    dry_bulb_temperature = 20. * units.degC
409
    wet_bulb_temperature = 18. * units.degC
410 View Code Duplication
    coeff = 6.1e-4 / units.kelvin
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
411
    psychrometric_rh = relative_humidity_wet_psychrometric(dry_bulb_temperature,
412
                                                           wet_bulb_temperature, p,
413
                                                           psychrometer_coefficient=coeff)
414
    assert_almost_equal(psychrometric_rh, 82.9701 * units.percent, 3)
415
416
417
def test_rh_mixing_ratio():
418
    """Tests relative humidity from mixing ratio."""
419
    p = 1013.25 * units.mbar
420
    temperature = 20. * units.degC
421
    w = 0.012
422
    rh = relative_humidity_from_mixing_ratio(w, temperature, p)
423
    assert_almost_equal(rh, 81.7219 * units.percent, 3)
424
425
426
def test_mixing_ratio_from_specific_humidity():
427
    """Tests mixing ratio from specific humidity."""
428
    q = 0.012
429
    w = mixing_ratio_from_specific_humidity(q)
430
    assert_almost_equal(w, 0.01215, 3)
431
432
433
def test_rh_specific_humidity():
434
    """Tests relative humidity from specific humidity."""
435
    p = 1013.25 * units.mbar
436
    temperature = 20. * units.degC
437 View Code Duplication
    q = 0.012
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
438
    rh = relative_humidity_from_specific_humidity(q, temperature, p)
439
    assert_almost_equal(rh, 82.7145 * units.percent, 3)
440
441
442
def test_cape_cin():
443
    """Tests the basic CAPE and CIN calculation."""
444
    p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
445
    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
446
    dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
447
    parcel_prof = parcel_profile(p, temperature[0], dewpoint[0])
448
    cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
449
    assert_almost_equal(cape, 58.0368212 * units('joule / kilogram'), 6)
450
    assert_almost_equal(cin, -89.8073512 * units('joule / kilogram'), 6)
451
452
453 View Code Duplication
def test_cape_cin_no_el():
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
454
    """Tests that CAPE works with no EL."""
455
    p = np.array([959., 779.2, 751.3, 724.3]) * units.mbar
456
    temperature = np.array([22.2, 14.6, 12., 9.4]) * units.celsius
457
    dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius
458
    parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC')
459
    cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
460
    assert_almost_equal(cape, 0.08750805 * units('joule / kilogram'), 6)
461
    assert_almost_equal(cin, -89.8073512 * units('joule / kilogram'), 6)
462
463
464
def test_cape_cin_no_lfc():
465
    """Tests that CAPE is zero with no LFC."""
466
    p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
467
    temperature = np.array([22.2, 24.6, 22., 20.4, 18., -10.]) * units.celsius
468 View Code Duplication
    dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
469
    parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC')
470
    cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
471
    assert_almost_equal(cape, 0.0 * units('joule / kilogram'), 6)
472
    assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 6)
473
474
475
def test_find_append_zero_crossings():
476
    """Tests finding and appending zero crossings of an x, y series."""
477
    x = np.arange(11) * units.hPa
478
    y = np.array([3, 2, 1, -1, 2, 2, 0, 1, 0, -1, 2]) * units.degC
479
    x2, y2 = _find_append_zero_crossings(x, y)
480
481
    x_truth = np.array([0., 1., 2., 2.5, 3., 3.33333333, 4., 5.,
482
                        6., 7., 8., 9., 9.33333333, 10.]) * units.hPa
483
    y_truth = np.array([3, 2, 1, 0, -1, 0, 2, 2, 0, 1, 0, -1, 0, 2]) * units.degC
484
    assert_array_almost_equal(x2, x_truth, 6)
485
    assert_almost_equal(y2, y_truth, 6)
486
487
488
def test_most_unstable_parcel():
489 View Code Duplication
    """Tests calculating the most unstable parcel."""
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
490
    levels = np.array([1000., 959., 867.9]) * units.mbar
491
    temperatures = np.array([18.2, 22.2, 17.4]) * units.celsius
492
    dewpoints = np.array([19., 19., 14.3]) * units.celsius
493
    ret = most_unstable_parcel(levels, temperatures, dewpoints, depth=100 * units.hPa)
494
    assert_almost_equal(ret[0], 959.0 * units.hPa, 6)
495
    assert_almost_equal(ret[1], 22.2 * units.degC, 6)
496
    assert_almost_equal(ret[2], 19.0 * units.degC, 6)
497
498
499
def test_isentropic_pressure():
500
    """Test calculation of isentropic pressure function."""
501
    lev = [100000., 95000., 90000., 85000.] * units.Pa
502
    tmp = np.ones((4, 5, 5))
503
    tmp[0, :] = 296.
504 View Code Duplication
    tmp[1, :] = 292.
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
505
    tmp[2, :] = 290
506
    tmp[3, :] = 288.
507
    tmpk = tmp * units.kelvin
508
    isentlev = [296.] * units.kelvin
509
    isentprs = isentropic_interpolation(isentlev, lev, tmpk)
510
    trueprs = 1000. * units.hPa
511
    assert_almost_equal(isentprs[0].shape, (1, 5, 5), 3)
512
    assert_almost_equal(isentprs[0], trueprs, 3)
513
514
515
def test_isentropic_pressure_p_increase():
516
    """Test calculation of isentropic pressure function, p increasing order."""
517
    lev = [85000, 90000., 95000., 100000.] * units.Pa
518
    tmp = np.ones((4, 5, 5))
519
    tmp[0, :] = 288.
520
    tmp[1, :] = 290.
521
    tmp[2, :] = 292.
522
    tmp[3, :] = 296.
523
    tmpk = tmp * units.kelvin
524
    isentlev = [296.] * units.kelvin
525 View Code Duplication
    isentprs = isentropic_interpolation(isentlev, lev, tmpk)
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
526
    trueprs = 1000. * units.hPa
527
    assert_almost_equal(isentprs[0], trueprs, 3)
528
529
530
def test_isentropic_pressure_adition_args():
531
    """Test calculation of isentropic pressure function, additional args."""
532
    lev = [100000., 95000., 90000., 85000.] * units.Pa
533
    tmp = np.ones((4, 5, 5))
534
    tmp[0, :] = 296.
535
    tmp[1, :] = 292.
536
    tmp[2, :] = 290.
537
    tmp[3, :] = 288.
538
    rh = np.ones((4, 5, 5))
539
    rh[0, :] = 100.
540 View Code Duplication
    rh[1, :] = 80.
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
541
    rh[2, :] = 40.
542
    rh[3, :] = 20.
543
    relh = rh * units.percent
544
    tmpk = tmp * units.kelvin
545
    isentlev = [296.] * units.kelvin
546
    isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh)
547
    truerh = 100. * units.percent
548
    assert_almost_equal(isentprs[1], truerh, 3)
549
550
551
def test_isentropic_pressure_tmp_out():
552
    """Test calculation of isentropic pressure function, temperature output."""
553
    lev = [100000., 95000., 90000., 85000.] * units.Pa
554
    tmp = np.ones((4, 5, 5))
555
    tmp[0, :] = 296.
556
    tmp[1, :] = 292.
557
    tmp[2, :] = 290.
558
    tmp[3, :] = 288.
559
    tmpk = tmp * units.kelvin
560
    isentlev = [296.] * units.kelvin
561 View Code Duplication
    isentprs = isentropic_interpolation(isentlev, lev, tmpk, tmpk_out=True)
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
562
    truetmp = 296. * units.kelvin
563
    assert_almost_equal(isentprs[1], truetmp, 3)
564
565
566
def test_isentropic_pressure_p_increase_rh_out():
567
    """Test calculation of isentropic pressure function, p increasing order."""
568
    lev = [85000., 90000., 95000., 100000.] * units.Pa
569
    tmp = np.ones((4, 5, 5))
570
    tmp[0, :] = 288.
571
    tmp[1, :] = 290.
572
    tmp[2, :] = 292.
573
    tmp[3, :] = 296.
574
    tmpk = tmp * units.kelvin
575
    rh = np.ones((4, 5, 5))
576
    rh[0, :] = 20.
577
    rh[1, :] = 40.
578
    rh[2, :] = 80.
579
    rh[3, :] = 100.
580
    relh = rh * units.percent
581
    isentlev = 296. * units.kelvin
582
    isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh)
583
    truerh = 100. * units.percent
584
    assert_almost_equal(isentprs[1], truerh, 3)
585
586
587
def test_isentropic_pressure_interp():
588
    """Test calculation of isentropic pressure function."""
589
    lev = [100000., 95000., 90000., 85000.] * units.Pa
590
    tmp = np.ones((4, 5, 5))
591
    tmp[0, :] = 296.
592
    tmp[1, :] = 292.
593
    tmp[2, :] = 290
594
    tmp[3, :] = 288.
595
    tmpk = tmp * units.kelvin
596
    isentlev = [296., 297] * units.kelvin
597
    isentprs = isentropic_interpolation(isentlev, lev, tmpk)
598
    trueprs = 936.18057 * units.hPa
599
    assert_almost_equal(isentprs[0][1], trueprs, 3)
600
601
602
def test_isentropic_pressure_adition_args_interp():
603
    """Test calculation of isentropic pressure function, additional args."""
604
    lev = [100000., 95000., 90000., 85000.] * units.Pa
605
    tmp = np.ones((4, 5, 5))
606
    tmp[0, :] = 296.
607
    tmp[1, :] = 292.
608
    tmp[2, :] = 290.
609
    tmp[3, :] = 288.
610
    rh = np.ones((4, 5, 5))
611
    rh[0, :] = 100.
612
    rh[1, :] = 80.
613
    rh[2, :] = 40.
614
    rh[3, :] = 20.
615
    relh = rh * units.percent
616
    tmpk = tmp * units.kelvin
617
    isentlev = [296., 297.] * units.kelvin
618
    isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh)
619
    truerh = 69.171 * units.percent
620
    assert_almost_equal(isentprs[1][1], truerh, 3)
621
622
623
def test_isentropic_pressure_tmp_out_interp():
624
    """Test calculation of isentropic pressure function, temperature output."""
625
    lev = [100000., 95000., 90000., 85000.] * units.Pa
626
    tmp = np.ones((4, 5, 5))
627
    tmp[0, :] = 296.
628
    tmp[1, :] = 292.
629
    tmp[2, :] = 290.
630
    tmp[3, :] = 288.
631
    tmpk = tmp * units.kelvin
632
    isentlev = [296., 297.] * units.kelvin
633
    isentprs = isentropic_interpolation(isentlev, lev, tmpk, tmpk_out=True)
634
    truetmp = 291.4579 * units.kelvin
635
    assert_almost_equal(isentprs[1][1], truetmp, 3)
636
637
638
def test_isentropic_pressure_data_bounds_error():
639
    """Test calculation of isentropic pressure function, error for data out of bounds."""
640
    lev = [100000., 95000., 90000., 85000.] * units.Pa
641
    tmp = np.ones((4, 5, 5))
642
    tmp[0, :] = 296.
643
    tmp[1, :] = 292.
644
    tmp[2, :] = 290.
645
    tmp[3, :] = 288.
646
    tmpk = tmp * units.kelvin
647
    isentlev = [296., 350.] * units.kelvin
648
    with pytest.raises(ValueError):
649
        isentropic_interpolation(isentlev, lev, tmpk)
650
651
652
def test_isentropic_pressure_4d():
653
    """Test calculation of isentropic pressure function."""
654
    lev = [100000., 95000., 90000., 85000.] * units.Pa
655
    tmp = np.ones((3, 4, 5, 5))
656
    tmp[:, 0, :] = 296.
657
    tmp[:, 1, :] = 292.
658
    tmp[:, 2, :] = 290
659
    tmp[:, 3, :] = 288.
660
    tmpk = tmp * units.kelvin
661
    rh = np.ones((3, 4, 5, 5))
662
    rh[:, 0, :] = 100.
663
    rh[:, 1, :] = 80.
664
    rh[:, 2, :] = 40.
665
    rh[:, 3, :] = 20.
666
    relh = rh * units.percent
667
    isentlev = [296., 297., 300.] * units.kelvin
668
    isentprs = isentropic_interpolation(isentlev, lev, tmpk, relh, axis=1)
669
    trueprs = 1000. * units.hPa
670
    trueprs2 = 936.18057 * units.hPa
671
    trueprs3 = 879.446 * units.hPa
672
    truerh = 69.171 * units.percent
673
    assert_almost_equal(isentprs[0].shape, (3, 3, 5, 5), 3)
674
    assert_almost_equal(isentprs[0][:, 0, :], trueprs, 3)
675
    assert_almost_equal(isentprs[0][:, 1, :], trueprs2, 3)
676
    assert_almost_equal(isentprs[0][:, 2, :], trueprs3, 3)
677
    assert_almost_equal(isentprs[1][:, 1, ], truerh, 3)
678
679
680
def test_surface_based_cape_cin():
681
    """Tests the surface-based CAPE and CIN calculation."""
682
    p = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
683
    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
684
    dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
685
    cape, cin = surface_based_cape_cin(p, temperature, dewpoint)
686
    assert_almost_equal(cape, 58.0368212 * units('joule / kilogram'), 6)
687
    assert_almost_equal(cin, -89.8073512 * units('joule / kilogram'), 6)
688
689
690
def test_most_unstable_cape_cin_surface():
691
    """Tests the most unstable CAPE/CIN calculation when surface is most unstable."""
692
    pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar
693
    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
694
    dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
695
    mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint)
696
    assert_almost_equal(mucape, 58.0368212 * units('joule / kilogram'), 6)
697
    assert_almost_equal(mucin, -89.8073512 * units('joule / kilogram'), 6)
698
699
700
def test_most_unstable_cape_cin():
701
    """Tests the most unstable CAPE/CIN calculation."""
702
    pressure = np.array([1000., 959., 867.9, 850., 825., 800.]) * units.mbar
703
    temperature = np.array([18.2, 22.2, 17.4, 10., 0., 15]) * units.celsius
704
    dewpoint = np.array([19., 19., 14.3, 0., -10., 0.]) * units.celsius
705
    mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint)
706
    assert_almost_equal(mucape, 157.07111 * units('joule / kilogram'), 4)
707
    assert_almost_equal(mucin, -15.74772 * units('joule / kilogram'), 4)
708
709
710
def test_mixed_parcel():
711
    """Tests the mixed parcel calculation."""
712
    pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa
713
    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC
714
    dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.degC
715
    parcel_pressure, parcel_temperature, parcel_dewpoint = mixed_parcel(pressure, temperature,
716
                                                                        dewpoint,
717
                                                                        depth=250 * units.hPa)
718
    assert_almost_equal(parcel_pressure, 959. * units.hPa, 6)
719
    assert_almost_equal(parcel_temperature, 28.7363771 * units.degC, 6)
720
    assert_almost_equal(parcel_dewpoint, 7.1534658 * units.degC, 6)
721
722
723
def test_mixed_layer():
724
    """Tests the mixed layer calculation."""
725
    pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.hPa
726
    temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.degC
727
    mixed_layer_temperature = mixed_layer(pressure, temperature, depth=250 * units.hPa)[0]
728
    assert_almost_equal(mixed_layer_temperature, 16.4024930 * units.degC, 6)
729