Completed
Push — master ( ab0a25...73ac07 )
by Ryan
16s
created

_check_radians()   B

Complexity

Conditions 3

Size

Total Lines 25

Duplication

Lines 0
Ratio 0 %

Importance

Changes 1
Bugs 0 Features 0
Metric Value
cc 3
c 1
b 0
f 0
dl 0
loc 25
rs 8.8571
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