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
|
|
|
"""Contains a collection of basic calculations. |
5
|
|
|
|
6
|
|
|
These include: |
7
|
|
|
|
8
|
|
|
* wind components |
9
|
|
|
* heat index |
10
|
|
|
* windchill |
11
|
|
|
""" |
12
|
|
|
|
13
|
|
|
from __future__ import division |
14
|
|
|
|
15
|
|
|
import warnings |
16
|
|
|
|
17
|
|
|
import numpy as np |
18
|
|
|
|
19
|
|
|
from ..constants import g, omega, Rd |
20
|
|
|
from ..package_tools import Exporter |
21
|
|
|
from ..units import atleast_1d, check_units, masked_array, units |
22
|
|
|
|
23
|
|
|
exporter = Exporter(globals()) |
24
|
|
|
|
25
|
|
|
|
26
|
|
|
@exporter.export |
27
|
|
|
def get_wind_speed(u, v): |
28
|
|
|
r"""Compute the wind speed from u and v-components. |
29
|
|
|
|
30
|
|
|
Parameters |
31
|
|
|
---------- |
32
|
|
|
u : array_like |
33
|
|
|
Wind component in the X (East-West) direction |
34
|
|
|
v : array_like |
35
|
|
|
Wind component in the Y (North-South) direction |
36
|
|
|
|
37
|
|
|
Returns |
38
|
|
|
------- |
39
|
|
|
wind speed: array_like |
40
|
|
|
The speed of the wind |
41
|
|
|
|
42
|
|
|
See Also |
43
|
|
|
-------- |
44
|
|
|
get_wind_components |
45
|
|
|
|
46
|
|
|
""" |
47
|
|
|
speed = np.sqrt(u * u + v * v) |
48
|
|
|
return speed |
49
|
|
|
|
50
|
|
|
|
51
|
|
|
@exporter.export |
52
|
|
|
def get_wind_dir(u, v): |
53
|
|
|
r"""Compute the wind direction from u and v-components. |
54
|
|
|
|
55
|
|
|
Parameters |
56
|
|
|
---------- |
57
|
|
|
u : array_like |
58
|
|
|
Wind component in the X (East-West) direction |
59
|
|
|
v : array_like |
60
|
|
|
Wind component in the Y (North-South) direction |
61
|
|
|
|
62
|
|
|
Returns |
63
|
|
|
------- |
64
|
|
|
wind direction: `pint.Quantity` |
65
|
|
|
The direction of the wind, specified as the direction from |
66
|
|
|
which it is blowing, with 0 being North. |
67
|
|
|
|
68
|
|
|
See Also |
69
|
|
|
-------- |
70
|
|
|
get_wind_components |
71
|
|
|
|
72
|
|
|
""" |
73
|
|
|
wdir = 90. * units.deg - np.arctan2(-v, -u) |
74
|
|
|
origshape = wdir.shape |
75
|
|
|
wdir = atleast_1d(wdir) |
76
|
|
|
wdir[wdir < 0] += 360. * units.deg |
77
|
|
|
return wdir.reshape(origshape) |
78
|
|
|
|
79
|
|
|
|
80
|
|
|
@exporter.export |
81
|
|
|
def get_wind_components(speed, wdir): |
82
|
|
|
r"""Calculate the U, V wind vector components from the speed and direction. |
83
|
|
|
|
84
|
|
|
Parameters |
85
|
|
|
---------- |
86
|
|
|
speed : array_like |
87
|
|
|
The wind speed (magnitude) |
88
|
|
|
wdir : array_like |
89
|
|
|
The wind direction, specified as the direction from which the wind is |
90
|
|
|
blowing, with 0 being North. |
91
|
|
|
|
92
|
|
|
Returns |
93
|
|
|
------- |
94
|
|
|
u, v : tuple of array_like |
95
|
|
|
The wind components in the X (East-West) and Y (North-South) |
96
|
|
|
directions, respectively. |
97
|
|
|
|
98
|
|
|
See Also |
99
|
|
|
-------- |
100
|
|
|
get_wind_speed |
101
|
|
|
get_wind_dir |
102
|
|
|
|
103
|
|
|
Examples |
104
|
|
|
-------- |
105
|
|
|
>>> from metpy.units import units |
106
|
|
|
>>> metpy.calc.get_wind_components(10. * units('m/s'), 225. * units.deg) |
107
|
|
|
(<Quantity(7.071067811865475, 'meter / second')>, |
108
|
|
|
<Quantity(7.071067811865477, 'meter / second')>) |
109
|
|
|
|
110
|
|
|
""" |
111
|
|
|
wdir = _check_radians(wdir, max_radians=4 * np.pi) |
112
|
|
|
u = -speed * np.sin(wdir) |
113
|
|
|
v = -speed * np.cos(wdir) |
114
|
|
|
return u, v |
115
|
|
|
|
116
|
|
|
|
117
|
|
|
@exporter.export |
118
|
|
|
@check_units(temperature='[temperature]', speed='[speed]') |
119
|
|
|
def windchill(temperature, speed, face_level_winds=False, mask_undefined=True): |
120
|
|
|
r"""Calculate the Wind Chill Temperature Index (WCTI). |
121
|
|
|
|
122
|
|
|
Calculates WCTI from the current temperature and wind speed using the formula |
123
|
|
|
outlined by the FCM [FCMR192003]_. |
124
|
|
|
|
125
|
|
|
Specifically, these formulas assume that wind speed is measured at |
126
|
|
|
10m. If, instead, the speeds are measured at face level, the winds |
127
|
|
|
need to be multiplied by a factor of 1.5 (this can be done by specifying |
128
|
|
|
`face_level_winds` as `True`.) |
129
|
|
|
|
130
|
|
|
Parameters |
131
|
|
|
---------- |
132
|
|
|
temperature : `pint.Quantity` |
133
|
|
|
The air temperature |
134
|
|
|
speed : `pint.Quantity` |
135
|
|
|
The wind speed at 10m. If instead the winds are at face level, |
136
|
|
|
`face_level_winds` should be set to `True` and the 1.5 multiplicative |
137
|
|
|
correction will be applied automatically. |
138
|
|
|
face_level_winds : bool, optional |
139
|
|
|
A flag indicating whether the wind speeds were measured at facial |
140
|
|
|
level instead of 10m, thus requiring a correction. Defaults to |
141
|
|
|
`False`. |
142
|
|
|
mask_undefined : bool, optional |
143
|
|
|
A flag indicating whether a masked array should be returned with |
144
|
|
|
values where wind chill is undefined masked. These are values where |
145
|
|
|
the temperature > 50F or wind speed <= 3 miles per hour. Defaults |
146
|
|
|
to `True`. |
147
|
|
|
|
148
|
|
|
Returns |
149
|
|
|
------- |
150
|
|
|
`pint.Quantity` |
151
|
|
|
The corresponding Wind Chill Temperature Index value(s) |
152
|
|
|
|
153
|
|
|
See Also |
154
|
|
|
-------- |
155
|
|
|
heat_index |
156
|
|
|
|
157
|
|
|
""" |
158
|
|
|
# Correct for lower height measurement of winds if necessary |
159
|
|
|
if face_level_winds: |
160
|
|
|
# No in-place so that we copy |
161
|
|
|
# noinspection PyAugmentAssignment |
162
|
|
|
speed = speed * 1.5 |
163
|
|
|
|
164
|
|
|
temp_limit, speed_limit = 10. * units.degC, 3 * units.mph |
165
|
|
|
speed_factor = speed.to('km/hr').magnitude ** 0.16 |
166
|
|
|
wcti = units.Quantity((0.6215 + 0.3965 * speed_factor) * temperature.to('degC').magnitude - |
167
|
|
|
11.37 * speed_factor + 13.12, units.degC).to(temperature.units) |
168
|
|
|
|
169
|
|
|
# See if we need to mask any undefined values |
170
|
|
|
if mask_undefined: |
171
|
|
|
mask = np.array((temperature > temp_limit) | (speed <= speed_limit)) |
172
|
|
|
if mask.any(): |
173
|
|
|
wcti = masked_array(wcti, mask=mask) |
174
|
|
|
|
175
|
|
|
return wcti |
176
|
|
|
|
177
|
|
|
|
178
|
|
|
@exporter.export |
179
|
|
|
@check_units('[temperature]') |
180
|
|
|
def heat_index(temperature, rh, mask_undefined=True): |
181
|
|
|
r"""Calculate the Heat Index from the current temperature and relative humidity. |
182
|
|
|
|
183
|
|
|
The implementation uses the formula outlined in [Rothfusz1990]_. This equation is a |
184
|
|
|
multi-variable least-squares regression of the values obtained in [Steadman1979]_. |
185
|
|
|
|
186
|
|
|
Parameters |
187
|
|
|
---------- |
188
|
|
|
temperature : `pint.Quantity` |
189
|
|
|
Air temperature |
190
|
|
|
rh : array_like |
191
|
|
|
The relative humidity expressed as a unitless ratio in the range [0, 1]. |
192
|
|
|
Can also pass a percentage if proper units are attached. |
193
|
|
|
|
194
|
|
|
Returns |
195
|
|
|
------- |
196
|
|
|
`pint.Quantity` |
197
|
|
|
The corresponding Heat Index value(s) |
198
|
|
|
|
199
|
|
|
Other Parameters |
200
|
|
|
---------------- |
201
|
|
|
mask_undefined : bool, optional |
202
|
|
|
A flag indicating whether a masked array should be returned with |
203
|
|
|
values where heat index is undefined masked. These are values where |
204
|
|
|
the temperature < 80F or relative humidity < 40 percent. Defaults |
205
|
|
|
to `True`. |
206
|
|
|
|
207
|
|
|
See Also |
208
|
|
|
-------- |
209
|
|
|
windchill |
210
|
|
|
|
211
|
|
|
""" |
212
|
|
|
delta = temperature - 0. * units.degF |
213
|
|
|
rh2 = rh * rh |
214
|
|
|
delta2 = delta * delta |
215
|
|
|
|
216
|
|
|
# Calculate the Heat Index -- constants converted for RH in [0, 1] |
217
|
|
|
hi = (-42.379 * units.degF + 2.04901523 * delta + |
218
|
|
|
1014.333127 * units.delta_degF * rh - 22.475541 * delta * rh - |
219
|
|
|
6.83783e-3 / units.delta_degF * delta2 - 5.481717e2 * units.delta_degF * rh2 + |
220
|
|
|
1.22874e-1 / units.delta_degF * delta2 * rh + 8.5282 * delta * rh2 - |
221
|
|
|
1.99e-2 / units.delta_degF * delta2 * rh2) |
222
|
|
|
|
223
|
|
|
# See if we need to mask any undefined values |
224
|
|
|
if mask_undefined: |
225
|
|
|
mask = np.array((temperature < 80. * units.degF) | (rh < 40 * units.percent)) |
226
|
|
|
if mask.any(): |
227
|
|
|
hi = masked_array(hi, mask=mask) |
228
|
|
|
|
229
|
|
|
return hi |
230
|
|
|
|
231
|
|
|
|
232
|
|
|
@exporter.export |
233
|
|
|
@check_units('[pressure]') |
234
|
|
|
def pressure_to_height_std(pressure): |
235
|
|
|
r"""Convert pressure data to heights using the U.S. standard atmosphere. |
236
|
|
|
|
237
|
|
|
The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61. |
238
|
|
|
|
239
|
|
|
Parameters |
240
|
|
|
---------- |
241
|
|
|
pressure : `pint.Quantity` |
242
|
|
|
Atmospheric pressure |
243
|
|
|
|
244
|
|
|
Returns |
245
|
|
|
------- |
246
|
|
|
`pint.Quantity` |
247
|
|
|
The corresponding height value(s) |
248
|
|
|
|
249
|
|
|
Notes |
250
|
|
|
----- |
251
|
|
|
.. math:: Z = \frac{T_0}{\Gamma}[1-\frac{p}{p_0}^\frac{R\Gamma}{g}] |
252
|
|
|
|
253
|
|
|
""" |
254
|
|
|
t0 = 288. * units.kelvin |
255
|
|
|
gamma = 6.5 * units('K/km') |
256
|
|
|
p0 = 1013.25 * units.mbar |
257
|
|
|
return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(Rd * gamma / g)) |
258
|
|
|
|
259
|
|
|
|
260
|
|
|
@exporter.export |
261
|
|
|
@check_units('[length]') |
262
|
|
|
def height_to_pressure_std(height): |
263
|
|
|
r"""Convert height data to pressures using the U.S. standard atmosphere. |
264
|
|
|
|
265
|
|
|
The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61. |
266
|
|
|
|
267
|
|
|
Parameters |
268
|
|
|
---------- |
269
|
|
|
height : `pint.Quantity` |
270
|
|
|
Atmospheric height |
271
|
|
|
|
272
|
|
|
Returns |
273
|
|
|
------- |
274
|
|
|
`pint.Quantity` |
275
|
|
|
The corresponding pressure value(s) |
276
|
|
|
|
277
|
|
|
Notes |
278
|
|
|
----- |
279
|
|
|
.. math:: p = p_0 e^{\frac{g}{R \Gamma} \text{ln}(1-\frac{Z \Gamma}{T_0})} |
280
|
|
|
|
281
|
|
|
""" |
282
|
|
|
t0 = 288. * units.kelvin |
283
|
|
|
gamma = 6.5 * units('K/km') |
284
|
|
|
p0 = 1013.25 * units.mbar |
285
|
|
|
return p0 * (1 - (gamma / t0) * height) ** (g / (Rd * gamma)) |
286
|
|
|
|
287
|
|
|
|
288
|
|
|
@exporter.export |
289
|
|
|
def coriolis_parameter(latitude): |
290
|
|
|
r"""Calculate the coriolis parameter at each point. |
291
|
|
|
|
292
|
|
|
The implementation uses the formula outlined in [Hobbs1977]_ pg.370-371. |
293
|
|
|
|
294
|
|
|
Parameters |
295
|
|
|
---------- |
296
|
|
|
latitude : array_like |
297
|
|
|
Latitude at each point |
298
|
|
|
|
299
|
|
|
Returns |
300
|
|
|
------- |
301
|
|
|
`pint.Quantity` |
302
|
|
|
The corresponding coriolis force at each point |
303
|
|
|
|
304
|
|
|
""" |
305
|
|
|
latitude = _check_radians(latitude, max_radians=np.pi / 2) |
306
|
|
|
return 2. * omega * np.sin(latitude) |
307
|
|
|
|
308
|
|
|
|
309
|
|
|
@exporter.export |
310
|
|
|
@check_units('[pressure]', '[length]') |
311
|
|
|
def add_height_to_pressure(pressure, height): |
312
|
|
|
r"""Calculate the pressure at a certain height above another pressure level. |
313
|
|
|
|
314
|
|
|
This assumes a standard atmosphere. |
315
|
|
|
|
316
|
|
|
Parameters |
317
|
|
|
---------- |
318
|
|
|
pressure : `pint.Quantity` |
319
|
|
|
Pressure level |
320
|
|
|
height : `pint.Quantity` |
321
|
|
|
Height above a pressure level |
322
|
|
|
|
323
|
|
|
Returns |
324
|
|
|
------- |
325
|
|
|
`pint.Quantity` |
326
|
|
|
The corresponding pressure value for the height above the pressure level |
327
|
|
|
|
328
|
|
|
See Also |
329
|
|
|
----- |
330
|
|
|
pressure_to_height_std, height_to_pressure_std, add_pressure_to_height |
331
|
|
|
|
332
|
|
|
""" |
333
|
|
|
pressure_level_height = pressure_to_height_std(pressure) |
334
|
|
|
return height_to_pressure_std(pressure_level_height + height) |
335
|
|
|
|
336
|
|
|
|
337
|
|
|
@exporter.export |
338
|
|
|
@check_units('[length]', '[pressure]') |
339
|
|
|
def add_pressure_to_height(height, pressure): |
340
|
|
|
r"""Calculate the height at a certain pressure above another height. |
341
|
|
|
|
342
|
|
|
This assumes a standard atmosphere. |
343
|
|
|
|
344
|
|
|
Parameters |
345
|
|
|
---------- |
346
|
|
|
height : `pint.Quantity` |
347
|
|
|
Height level |
348
|
|
|
pressure : `pint.Quantity` |
349
|
|
|
Pressure above height level |
350
|
|
|
|
351
|
|
|
Returns |
352
|
|
|
------- |
353
|
|
|
`pint.Quantity` |
354
|
|
|
The corresponding height value for the pressure above the height level |
355
|
|
|
|
356
|
|
|
See Also |
357
|
|
|
----- |
358
|
|
|
pressure_to_height_std, height_to_pressure_std, add_height_to_pressure |
359
|
|
|
|
360
|
|
|
""" |
361
|
|
|
pressure_at_height = height_to_pressure_std(height) |
362
|
|
|
return pressure_to_height_std(pressure_at_height - pressure) |
363
|
|
|
|
364
|
|
|
|
365
|
|
|
@exporter.export |
366
|
|
|
@check_units('[dimensionless]', '[pressure]', '[pressure]') |
367
|
|
|
def sigma_to_pressure(sigma, psfc, ptop): |
368
|
|
|
r"""Calculate pressure from sigma values. |
369
|
|
|
|
370
|
|
|
Parameters |
371
|
|
|
---------- |
372
|
|
|
sigma : ndarray |
373
|
|
|
The sigma levels to be converted to pressure levels. |
374
|
|
|
|
375
|
|
|
psfc : `pint.Quantity` |
376
|
|
|
The surface pressure value. |
377
|
|
|
|
378
|
|
|
ptop : `pint.Quantity` |
379
|
|
|
The pressure value at the top of the model domain. |
380
|
|
|
|
381
|
|
|
Returns |
382
|
|
|
------- |
383
|
|
|
`pint.Quantity` |
384
|
|
|
The pressure values at the given sigma levels. |
385
|
|
|
|
386
|
|
|
Notes |
387
|
|
|
----- |
388
|
|
|
Sigma definition adapted from [Philips1957]_. |
389
|
|
|
|
390
|
|
|
.. math:: p = \sigma * (p_{sfc} - p_{top}) + p_{top} |
391
|
|
|
|
392
|
|
|
* :math:`p` is pressure at a given `\sigma` level |
393
|
|
|
* :math:`\sigma` is non-dimensional, scaled pressure |
394
|
|
|
* :math:`p_{sfc}` is pressure at the surface or model floor |
395
|
|
|
* :math:`p_{top}` is pressure at the top of the model domain |
396
|
|
|
|
397
|
|
|
""" |
398
|
|
|
if np.any(sigma < 0) or np.any(sigma > 1): |
399
|
|
|
raise ValueError('Sigma values should be bounded by 0 and 1') |
400
|
|
|
|
401
|
|
|
if psfc.magnitude < 0 or ptop.magnitude < 0: |
402
|
|
|
raise ValueError('Pressure input should be non-negative') |
403
|
|
|
|
404
|
|
|
return sigma * (psfc - ptop) + ptop |
405
|
|
|
|
406
|
|
|
|
407
|
|
|
def _check_radians(value, max_radians=2 * np.pi): |
408
|
|
|
"""Input validation of values that could be in degrees instead of radians. |
409
|
|
|
|
410
|
|
|
Parameters |
411
|
|
|
---------- |
412
|
|
|
value : `pint.Quantity` |
413
|
|
|
The input value to check. |
414
|
|
|
|
415
|
|
|
max_radians : float |
416
|
|
|
Maximum absolute value of radians before warning. |
417
|
|
|
|
418
|
|
|
Returns |
419
|
|
|
------- |
420
|
|
|
`pint.Quantity` |
421
|
|
|
The input value |
422
|
|
|
|
423
|
|
|
""" |
424
|
|
|
try: |
425
|
|
|
value = value.to('radians').m |
426
|
|
|
except AttributeError: |
427
|
|
|
pass |
428
|
|
|
if np.greater(np.nanmax(np.abs(value)), max_radians): |
429
|
|
|
warnings.warn('Input over {} radians. ' |
430
|
|
|
'Ensure proper units are given.'.format(max_radians)) |
431
|
|
|
return value |
432
|
|
|
|