1# Copyright (c) 2008,2015,2016,2017,2018,2019 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
6These include:
7
8* wind components
9* heat index
10* windchill
11"""
12import contextlib
13from itertools import product
14import warnings
15
16import numpy as np
17from scipy.ndimage import gaussian_filter
18
19from .. import constants as mpconsts
20from ..package_tools import Exporter
21from ..units import check_units, masked_array, units
22from ..xarray import preprocess_and_wrap
23
24exporter = Exporter(globals())
25
26# The following variables are constants for a standard atmosphere
27t0 = units.Quantity(288., 'kelvin')
28p0 = units.Quantity(1013.25, 'hPa')
29gamma = units.Quantity(6.5, 'K/km')
30
31
32@exporter.export
33@preprocess_and_wrap(wrap_like='u')
34@check_units('[speed]', '[speed]')
35def wind_speed(u, v):
36    r"""Compute the wind speed from u and v-components.
37
38    Parameters
39    ----------
40    u : `pint.Quantity`
41        Wind component in the X (East-West) direction
42    v : `pint.Quantity`
43        Wind component in the Y (North-South) direction
44
45    Returns
46    -------
47    wind speed: `pint.Quantity`
48        Speed of the wind
49
50    See Also
51    --------
52    wind_components
53
54    """
55    speed = np.sqrt(u * u + v * v)
56    return speed
57
58
59@exporter.export
60@preprocess_and_wrap(wrap_like='u')
61@check_units('[speed]', '[speed]')
62def wind_direction(u, v, convention='from'):
63    r"""Compute the wind direction from u and v-components.
64
65    Parameters
66    ----------
67    u : `pint.Quantity`
68        Wind component in the X (East-West) direction
69    v : `pint.Quantity`
70        Wind component in the Y (North-South) direction
71    convention : str
72        Convention to return direction; 'from' returns the direction the wind is coming from
73        (meteorological convention), 'to' returns the direction the wind is going towards
74        (oceanographic convention), default is 'from'.
75
76    Returns
77    -------
78    direction: `pint.Quantity`
79        The direction of the wind in intervals [0, 360] degrees, with 360 being North,
80        direction defined by the convention kwarg.
81
82    See Also
83    --------
84    wind_components
85
86    Notes
87    -----
88    In the case of calm winds (where `u` and `v` are zero), this function returns a direction
89    of 0.
90
91    """
92    wdir = units.Quantity(90., 'deg') - np.arctan2(-v, -u)
93    origshape = wdir.shape
94    wdir = np.atleast_1d(wdir)
95
96    # Handle oceanographic convection
97    if convention == 'to':
98        wdir -= units.Quantity(180., 'deg')
99    elif convention not in ('to', 'from'):
100        raise ValueError('Invalid kwarg for "convention". Valid options are "from" or "to".')
101
102    mask = wdir <= 0
103    if np.any(mask):
104        wdir[mask] += units.Quantity(360., 'deg')
105    # avoid unintended modification of `pint.Quantity` by direct use of magnitude
106    calm_mask = (np.asarray(u.magnitude) == 0.) & (np.asarray(v.magnitude) == 0.)
107    # np.any check required for legacy numpy which treats 0-d False boolean index as zero
108    if np.any(calm_mask):
109        wdir[calm_mask] = units.Quantity(0., 'deg')
110    return wdir.reshape(origshape).to('degrees')
111
112
113@exporter.export
114@preprocess_and_wrap(wrap_like=('speed', 'speed'))
115@check_units('[speed]')
116def wind_components(speed, wind_direction):
117    r"""Calculate the U, V wind vector components from the speed and direction.
118
119    Parameters
120    ----------
121    speed : `pint.Quantity`
122        Wind speed (magnitude)
123    wind_direction : `pint.Quantity`
124        Wind direction, specified as the direction from which the wind is
125        blowing (0-2 pi radians or 0-360 degrees), with 360 degrees being North.
126
127    Returns
128    -------
129    u, v : tuple of `pint.Quantity`
130        The wind components in the X (East-West) and Y (North-South)
131        directions, respectively.
132
133    See Also
134    --------
135    wind_speed
136    wind_direction
137
138    Examples
139    --------
140    >>> from metpy.units import units
141    >>> metpy.calc.wind_components(10. * units('m/s'), 225. * units.deg)
142     (<Quantity(7.07106781, 'meter / second')>, <Quantity(7.07106781, 'meter / second')>)
143
144    .. versionchanged:: 1.0
145       Renamed ``wdir`` parameter to ``wind_direction``
146
147    """
148    wind_direction = _check_radians(wind_direction, max_radians=4 * np.pi)
149    u = -speed * np.sin(wind_direction)
150    v = -speed * np.cos(wind_direction)
151    return u, v
152
153
154@exporter.export
155@preprocess_and_wrap(wrap_like='temperature')
156@check_units(temperature='[temperature]', speed='[speed]')
157def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
158    r"""Calculate the Wind Chill Temperature Index (WCTI).
159
160    Calculates WCTI from the current temperature and wind speed using the formula
161    outlined by the FCM [FCMR192003]_.
162
163    Specifically, these formulas assume that wind speed is measured at
164    10m.  If, instead, the speeds are measured at face level, the winds
165    need to be multiplied by a factor of 1.5 (this can be done by specifying
166    `face_level_winds` as `True`).
167
168    Parameters
169    ----------
170    temperature : `pint.Quantity`
171        Air temperature
172    speed : `pint.Quantity`
173        Wind speed at 10m. If instead the winds are at face level,
174        `face_level_winds` should be set to `True` and the 1.5 multiplicative
175        correction will be applied automatically.
176    face_level_winds : bool, optional
177        A flag indicating whether the wind speeds were measured at facial
178        level instead of 10m, thus requiring a correction.  Defaults to
179        `False`.
180    mask_undefined : bool, optional
181        A flag indicating whether a masked array should be returned with
182        values where wind chill is undefined masked.  These are values where
183        the temperature > 50F or wind speed <= 3 miles per hour. Defaults
184        to `True`.
185
186    Returns
187    -------
188    `pint.Quantity`
189        Corresponding Wind Chill Temperature Index value(s)
190
191    See Also
192    --------
193    heat_index, apparent_temperature
194
195    """
196    # Correct for lower height measurement of winds if necessary
197    if face_level_winds:
198        # No in-place so that we copy
199        # noinspection PyAugmentAssignment
200        speed = speed * 1.5
201
202    temp_limit, speed_limit = units.Quantity(10., 'degC'), units.Quantity(3, 'mph')
203    speed_factor = speed.to('km/hr').magnitude ** 0.16
204    wcti = units.Quantity((0.6215 + 0.3965 * speed_factor) * temperature.to('degC').magnitude
205                          - 11.37 * speed_factor + 13.12, units.degC).to(temperature.units)
206
207    # See if we need to mask any undefined values
208    if mask_undefined:
209        mask = np.array((temperature > temp_limit) | (speed <= speed_limit))
210        if mask.any():
211            wcti = masked_array(wcti, mask=mask)
212
213    return wcti
214
215
216@exporter.export
217@preprocess_and_wrap(wrap_like='temperature')
218@check_units('[temperature]')
219def heat_index(temperature, relative_humidity, mask_undefined=True):
220    r"""Calculate the Heat Index from the current temperature and relative humidity.
221
222    The implementation uses the formula outlined in [Rothfusz1990]_, which is a
223    multi-variable least-squares regression of the values obtained in [Steadman1979]_.
224    Additional conditional corrections are applied to match what the National
225    Weather Service operationally uses. See Figure 3 of [Anderson2013]_ for a
226    depiction of this algorithm and further discussion.
227
228    Parameters
229    ----------
230    temperature : `pint.Quantity`
231        Air temperature
232    relative_humidity : `pint.Quantity`
233        The relative humidity expressed as a unitless ratio in the range [0, 1].
234        Can also pass a percentage if proper units are attached.
235
236    Returns
237    -------
238    `pint.Quantity`
239        Corresponding Heat Index value(s)
240
241    Other Parameters
242    ----------------
243    mask_undefined : bool, optional
244        A flag indicating whether a masked array should be returned with
245        values masked where the temperature < 80F. Defaults to `True`.
246
247
248    .. versionchanged:: 1.0
249       Renamed ``rh`` parameter to ``relative_humidity``
250
251    See Also
252    --------
253    windchill, apparent_temperature
254
255    """
256    temperature = np.atleast_1d(temperature)
257    relative_humidity = np.atleast_1d(relative_humidity)
258    # assign units to relative_humidity if they currently are not present
259    if not hasattr(relative_humidity, 'units'):
260        relative_humidity = units.Quantity(relative_humidity, 'dimensionless')
261    delta = temperature.to(units.degF) - units.Quantity(0., 'degF')
262    rh2 = relative_humidity * relative_humidity
263    delta2 = delta * delta
264
265    # Simplifed Heat Index -- constants converted for relative_humidity in [0, 1]
266    a = (units.Quantity(-10.3, 'degF') + 1.1 * delta
267         + units.Quantity(4.7, 'delta_degF') * relative_humidity)
268
269    # More refined Heat Index -- constants converted for relative_humidity in [0, 1]
270    b = (units.Quantity(-42.379, 'degF')
271         + 2.04901523 * delta
272         + units.Quantity(1014.333127, 'delta_degF') * relative_humidity
273         - 22.475541 * delta * relative_humidity
274         - units.Quantity(6.83783e-3, '1/delta_degF') * delta2
275         - units.Quantity(5.481717e2, 'delta_degF') * rh2
276         + units.Quantity(1.22874e-1, '1/delta_degF') * delta2 * relative_humidity
277         + 8.5282 * delta * rh2
278         - units.Quantity(1.99e-2, '1/delta_degF') * delta2 * rh2)
279
280    # Create return heat index
281    hi = units.Quantity(np.full(np.shape(temperature), np.nan), 'degF')
282    # Retain masked status of temperature with resulting heat index
283    if hasattr(temperature, 'mask'):
284        hi = masked_array(hi)
285
286    # If T <= 40F, Heat Index is T
287    sel = (temperature <= units.Quantity(40., 'degF'))
288    if np.any(sel):
289        hi[sel] = temperature[sel].to(units.degF)
290
291    # If a < 79F and hi is unset, Heat Index is a
292    sel = (a < units.Quantity(79., 'degF')) & np.isnan(hi)
293    if np.any(sel):
294        hi[sel] = a[sel]
295
296    # Use b now for anywhere hi has yet to be set
297    sel = np.isnan(hi)
298    if np.any(sel):
299        hi[sel] = b[sel]
300
301    # Adjustment for RH <= 13% and 80F <= T <= 112F
302    sel = ((relative_humidity <= units.Quantity(13., 'percent'))
303           & (temperature >= units.Quantity(80., 'degF'))
304           & (temperature <= units.Quantity(112., 'degF')))
305    if np.any(sel):
306        rh15adj = ((13. - relative_humidity[sel] * 100.) / 4.
307                   * np.sqrt((units.Quantity(17., 'delta_degF')
308                              - np.abs(delta[sel] - units.Quantity(95., 'delta_degF')))
309                             / units.Quantity(17., '1/delta_degF')))
310        hi[sel] = hi[sel] - rh15adj
311
312    # Adjustment for RH > 85% and 80F <= T <= 87F
313    sel = ((relative_humidity > units.Quantity(85., 'percent'))
314           & (temperature >= units.Quantity(80., 'degF'))
315           & (temperature <= units.Quantity(87., 'degF')))
316    if np.any(sel):
317        rh85adj = (0.02 * (relative_humidity[sel] * 100. - 85.)
318                   * (units.Quantity(87., 'delta_degF') - delta[sel]))
319        hi[sel] = hi[sel] + rh85adj
320
321    # See if we need to mask any undefined values
322    if mask_undefined:
323        mask = np.array(temperature < units.Quantity(80., 'degF'))
324        if mask.any():
325            hi = masked_array(hi, mask=mask)
326    return hi
327
328
329@exporter.export
330@preprocess_and_wrap(wrap_like='temperature')
331@check_units(temperature='[temperature]', speed='[speed]')
332def apparent_temperature(temperature, relative_humidity, speed, face_level_winds=False,
333                         mask_undefined=True):
334    r"""Calculate the current apparent temperature.
335
336    Calculates the current apparent temperature based on the wind chill or heat index
337    as appropriate for the current conditions. Follows [NWS10201]_.
338
339    Parameters
340    ----------
341    temperature : `pint.Quantity`
342        Air temperature
343    relative_humidity : `pint.Quantity`
344        Relative humidity expressed as a unitless ratio in the range [0, 1].
345        Can also pass a percentage if proper units are attached.
346    speed : `pint.Quantity`
347        Wind speed at 10m.  If instead the winds are at face level,
348        `face_level_winds` should be set to `True` and the 1.5 multiplicative
349        correction will be applied automatically.
350    face_level_winds : bool, optional
351        A flag indicating whether the wind speeds were measured at facial
352        level instead of 10m, thus requiring a correction.  Defaults to
353        `False`.
354    mask_undefined : bool, optional
355        A flag indicating whether a masked array should be returned with
356        values where wind chill or heat_index is undefined masked. For wind
357        chill, these are values where the temperature > 50F or
358        wind speed <= 3 miles per hour. For heat index, these are values
359        where the temperature < 80F.
360        Defaults to `True`.
361
362    Returns
363    -------
364    `pint.Quantity`
365        Corresponding apparent temperature value(s)
366
367
368    .. versionchanged:: 1.0
369       Renamed ``rh`` parameter to ``relative_humidity``
370
371    See Also
372    --------
373    heat_index, windchill
374
375    """
376    is_not_scalar = isinstance(temperature.m, (list, tuple, np.ndarray))
377
378    temperature = np.atleast_1d(temperature)
379    relative_humidity = np.atleast_1d(relative_humidity)
380    speed = np.atleast_1d(speed)
381
382    # NB: mask_defined=True is needed to know where computed values exist
383    wind_chill_temperature = windchill(temperature, speed, face_level_winds=face_level_winds,
384                                       mask_undefined=True).to(temperature.units)
385
386    heat_index_temperature = heat_index(temperature, relative_humidity,
387                                        mask_undefined=True).to(temperature.units)
388
389    # Combine the heat index and wind chill arrays (no point has a value in both)
390    # NB: older numpy.ma.where does not return a masked array
391    app_temperature = masked_array(
392        np.ma.where(masked_array(wind_chill_temperature).mask,
393                    heat_index_temperature.to(temperature.units),
394                    wind_chill_temperature.to(temperature.units)
395                    ), temperature.units)
396
397    # If mask_undefined is False, then set any masked values to the temperature
398    if not mask_undefined:
399        app_temperature[app_temperature.mask] = temperature[app_temperature.mask]
400
401    # If no values are masked and provided temperature does not have a mask
402    # we should return a non-masked array
403    if not np.any(app_temperature.mask) and not hasattr(temperature, 'mask'):
404        app_temperature = units.Quantity(np.array(app_temperature.m), temperature.units)
405
406    if is_not_scalar:
407        return app_temperature
408    else:
409        return np.atleast_1d(app_temperature)[0]
410
411
412@exporter.export
413@preprocess_and_wrap(wrap_like='pressure')
414@check_units('[pressure]')
415def pressure_to_height_std(pressure):
416    r"""Convert pressure data to height using the U.S. standard atmosphere [NOAA1976]_.
417
418    The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61.
419
420    Parameters
421    ----------
422    pressure : `pint.Quantity`
423        Atmospheric pressure
424
425    Returns
426    -------
427    `pint.Quantity`
428        Corresponding height value(s)
429
430    Notes
431    -----
432    .. math:: Z = \frac{T_0}{\Gamma}[1-\frac{p}{p_0}^\frac{R\Gamma}{g}]
433
434    """
435    return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(
436        mpconsts.Rd * gamma / mpconsts.g))
437
438
439@exporter.export
440@preprocess_and_wrap(wrap_like='height')
441@check_units('[length]')
442def height_to_geopotential(height):
443    r"""Compute geopotential for a given height above sea level.
444
445    Calculates the geopotential from height above mean sea level using the following formula,
446    which is derived from the definition of geopotential as given in [Hobbs2006]_ Pg. 69 Eq
447    3.21, along with an approximation for variation of gravity with altitude:
448
449    .. math:: \Phi = \frac{g R_e z}{R_e + z}
450
451    (where :math:`\Phi` is geopotential, :math:`z` is height, :math:`R_e` is average Earth
452    radius, and :math:`g` is standard gravity).
453
454
455    Parameters
456    ----------
457    height : `pint.Quantity`
458        Height above sea level
459
460    Returns
461    -------
462    `pint.Quantity`
463        Corresponding geopotential value(s)
464
465    Examples
466    --------
467    >>> import metpy.calc
468    >>> from metpy.units import units
469    >>> height = np.linspace(0, 10000, num=11) * units.m
470    >>> geopot = metpy.calc.height_to_geopotential(height)
471    >>> geopot
472    <Quantity([     0.           9805.11097983 19607.1448853  29406.10316465
473    39201.98726524 48994.79863351 58784.53871501 68571.20895435
474    78354.81079527 88135.34568058 97912.81505219], 'meter ** 2 / second ** 2')>
475
476    Notes
477    -----
478    This calculation approximates :math:`g(z)` as
479
480    .. math:: g(z) = g_0 \left( \frac{R_e}{R_e + z} \right)^2
481
482    where :math:`g_0` is standard gravity. It thereby accounts for the average effects of
483    centrifugal force on apparent gravity, but neglects latitudinal variations due to
484    centrifugal force and Earth's eccentricity.
485
486    (Prior to MetPy v0.11, this formula instead calculated :math:`g(z)` from Newton's Law of
487    Gravitation assuming a spherical Earth and no centrifugal force effects).
488
489    See Also
490    --------
491    geopotential_to_height
492
493    """
494    return (mpconsts.g * mpconsts.Re * height) / (mpconsts.Re + height)
495
496
497@exporter.export
498@preprocess_and_wrap(wrap_like='geopotential')
499@check_units('[length] ** 2 / [time] ** 2')
500def geopotential_to_height(geopotential):
501    r"""Compute height above sea level from a given geopotential.
502
503    Calculates the height above mean sea level from geopotential using the following formula,
504    which is derived from the definition of geopotential as given in [Hobbs2006]_ Pg. 69 Eq
505    3.21, along with an approximation for variation of gravity with altitude:
506
507    .. math:: z = \frac{\Phi R_e}{gR_e - \Phi}
508
509    (where :math:`\Phi` is geopotential, :math:`z` is height, :math:`R_e` is average Earth
510    radius, and :math:`g` is standard gravity).
511
512
513    Parameters
514    ----------
515    geopotential : `pint.Quantity`
516        Geopotential
517
518    Returns
519    -------
520    `pint.Quantity`
521        Corresponding value(s) of height above sea level
522
523    Examples
524    --------
525    >>> import metpy.calc
526    >>> from metpy.units import units
527    >>> height = np.linspace(0, 10000, num=11) * units.m
528    >>> geopot = metpy.calc.height_to_geopotential(height)
529    >>> geopot
530    <Quantity([     0.           9805.11097983 19607.1448853  29406.10316465
531    39201.98726524 48994.79863351 58784.53871501 68571.20895435
532    78354.81079527 88135.34568058 97912.81505219], 'meter ** 2 / second ** 2')>
533    >>> height = metpy.calc.geopotential_to_height(geopot)
534    >>> height
535    <Quantity([     0.   1000.   2000.   3000.   4000.   5000.   6000.   7000.   8000.
536    9000.  10000.], 'meter')>
537
538    Notes
539    -----
540    This calculation approximates :math:`g(z)` as
541
542    .. math:: g(z) = g_0 \left( \frac{R_e}{R_e + z} \right)^2
543
544    where :math:`g_0` is standard gravity. It thereby accounts for the average effects of
545    centrifugal force on apparent gravity, but neglects latitudinal variations due to
546    centrifugal force and Earth's eccentricity.
547
548    (Prior to MetPy v0.11, this formula instead calculated :math:`g(z)` from Newton's Law of
549    Gravitation assuming a spherical Earth and no centrifugal force effects.)
550
551    .. versionchanged:: 1.0
552       Renamed ``geopot`` parameter to ``geopotential``
553
554    See Also
555    --------
556    height_to_geopotential
557
558    """
559    return (geopotential * mpconsts.Re) / (mpconsts.g * mpconsts.Re - geopotential)
560
561
562@exporter.export
563@preprocess_and_wrap(wrap_like='height')
564@check_units('[length]')
565def height_to_pressure_std(height):
566    r"""Convert height data to pressures using the U.S. standard atmosphere [NOAA1976]_.
567
568    The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61.
569
570    Parameters
571    ----------
572    height : `pint.Quantity`
573        Atmospheric height
574
575    Returns
576    -------
577    `pint.Quantity`
578        Corresponding pressure value(s)
579
580    Notes
581    -----
582    .. math:: p = p_0 e^{\frac{g}{R \Gamma} \text{ln}(1-\frac{Z \Gamma}{T_0})}
583
584    """
585    return p0 * (1 - (gamma / t0) * height) ** (mpconsts.g / (mpconsts.Rd * gamma))
586
587
588@exporter.export
589@preprocess_and_wrap(wrap_like='latitude')
590def coriolis_parameter(latitude):
591    r"""Calculate the coriolis parameter at each point.
592
593    The implementation uses the formula outlined in [Hobbs1977]_ pg.370-371.
594
595    Parameters
596    ----------
597    latitude : array_like
598        Latitude at each point
599
600    Returns
601    -------
602    `pint.Quantity`
603        Corresponding coriolis force at each point
604
605    """
606    latitude = _check_radians(latitude, max_radians=np.pi / 2)
607    return (2. * mpconsts.omega * np.sin(latitude)).to('1/s')
608
609
610@exporter.export
611@preprocess_and_wrap(wrap_like='pressure')
612@check_units('[pressure]', '[length]')
613def add_height_to_pressure(pressure, height):
614    r"""Calculate the pressure at a certain height above another pressure level.
615
616    This assumes a standard atmosphere [NOAA1976]_.
617
618    Parameters
619    ----------
620    pressure : `pint.Quantity`
621        Pressure level
622    height : `pint.Quantity`
623        Height above a pressure level
624
625    Returns
626    -------
627    `pint.Quantity`
628        Corresponding pressure value for the height above the pressure level
629
630    See Also
631    --------
632    pressure_to_height_std, height_to_pressure_std, add_pressure_to_height
633
634    """
635    pressure_level_height = pressure_to_height_std(pressure)
636    return height_to_pressure_std(pressure_level_height + height)
637
638
639@exporter.export
640@preprocess_and_wrap(wrap_like='height')
641@check_units('[length]', '[pressure]')
642def add_pressure_to_height(height, pressure):
643    r"""Calculate the height at a certain pressure above another height.
644
645    This assumes a standard atmosphere [NOAA1976]_.
646
647    Parameters
648    ----------
649    height : `pint.Quantity`
650        Height level
651    pressure : `pint.Quantity`
652        Pressure above height level
653
654    Returns
655    -------
656    `pint.Quantity`
657        The corresponding height value for the pressure above the height level
658
659    See Also
660    --------
661    pressure_to_height_std, height_to_pressure_std, add_height_to_pressure
662
663    """
664    pressure_at_height = height_to_pressure_std(height)
665    return pressure_to_height_std(pressure_at_height - pressure)
666
667
668@exporter.export
669@preprocess_and_wrap(wrap_like='sigma')
670@check_units('[dimensionless]', '[pressure]', '[pressure]')
671def sigma_to_pressure(sigma, pressure_sfc, pressure_top):
672    r"""Calculate pressure from sigma values.
673
674    Parameters
675    ----------
676    sigma : ndarray
677        Sigma levels to be converted to pressure levels
678
679    pressure_sfc : `pint.Quantity`
680        Surface pressure value
681
682    pressure_top : `pint.Quantity`
683        Pressure value at the top of the model domain
684
685    Returns
686    -------
687    `pint.Quantity`
688        Pressure values at the given sigma levels
689
690    Notes
691    -----
692    Sigma definition adapted from [Philips1957]_:
693
694    .. math:: p = \sigma * (p_{sfc} - p_{top}) + p_{top}
695
696    * :math:`p` is pressure at a given `\sigma` level
697    * :math:`\sigma` is non-dimensional, scaled pressure
698    * :math:`p_{sfc}` is pressure at the surface or model floor
699    * :math:`p_{top}` is pressure at the top of the model domain
700
701    .. versionchanged:: 1.0
702       Renamed ``psfc``, ``ptop`` parameters to ``pressure_sfc``, ``pressure_top``
703
704    """
705    if np.any(sigma < 0) or np.any(sigma > 1):
706        raise ValueError('Sigma values should be bounded by 0 and 1')
707
708    if pressure_sfc.magnitude < 0 or pressure_top.magnitude < 0:
709        raise ValueError('Pressure input should be non-negative')
710
711    return sigma * (pressure_sfc - pressure_top) + pressure_top
712
713
714@exporter.export
715@preprocess_and_wrap(wrap_like='scalar_grid', match_unit=True, to_magnitude=True)
716def smooth_gaussian(scalar_grid, n):
717    """Filter with normal distribution of weights.
718
719    Parameters
720    ----------
721    scalar_grid : `pint.Quantity`
722        Some n-dimensional scalar grid. If more than two axes, smoothing
723        is only done across the last two.
724
725    n : int
726        Degree of filtering
727
728    Returns
729    -------
730    `pint.Quantity`
731        The filtered 2D scalar grid
732
733    Notes
734    -----
735    This function is a close replication of the GEMPAK function ``GWFS``,
736    but is not identical.  The following notes are incorporated from
737    the GEMPAK source code:
738
739    This function smoothes a scalar grid using a moving average
740    low-pass filter whose weights are determined by the normal
741    (Gaussian) probability distribution function for two dimensions.
742    The weight given to any grid point within the area covered by the
743    moving average for a target grid point is proportional to:
744
745    .. math:: e^{-D^2}
746
747    where D is the distance from that point to the target point divided
748    by the standard deviation of the normal distribution.  The value of
749    the standard deviation is determined by the degree of filtering
750    requested.  The degree of filtering is specified by an integer.
751    This integer is the number of grid increments from crest to crest
752    of the wave for which the theoretical response is 1/e = .3679.  If
753    the grid increment is called delta_x, and the value of this integer
754    is represented by N, then the theoretical filter response function
755    value for the N * delta_x wave will be 1/e.  The actual response
756    function will be greater than the theoretical value.
757
758    The larger N is, the more severe the filtering will be, because the
759    response function for all wavelengths shorter than N * delta_x
760    will be less than 1/e.  Furthermore, as N is increased, the slope
761    of the filter response function becomes more shallow; so, the
762    response at all wavelengths decreases, but the amount of decrease
763    lessens with increasing wavelength.  (The theoretical response
764    function can be obtained easily--it is the Fourier transform of the
765    weight function described above.)
766
767    The area of the patch covered by the moving average varies with N.
768    As N gets bigger, the smoothing gets stronger, and weight values
769    farther from the target grid point are larger because the standard
770    deviation of the normal distribution is bigger.  Thus, increasing
771    N has the effect of expanding the moving average window as well as
772    changing the values of weights.  The patch is a square covering all
773    points whose weight values are within two standard deviations of the
774    mean of the two dimensional normal distribution.
775
776    The key difference between GEMPAK's GWFS and this function is that,
777    in GEMPAK, the leftover weight values representing the fringe of the
778    distribution are applied to the target grid point.  In this
779    function, the leftover weights are not used.
780
781    When this function is invoked, the first argument is the grid to be
782    smoothed, the second is the value of N as described above:
783
784                        GWFS ( S, N )
785
786    where N > 1.  If N <= 1, N = 2 is assumed.  For example, if N = 4,
787    then the 4 delta x wave length is passed with approximate response
788    1/e.
789
790    """
791    # Compute standard deviation in a manner consistent with GEMPAK
792    n = int(round(n))
793    if n < 2:
794        n = 2
795    sgma = n / (2 * np.pi)
796
797    # Construct sigma sequence so smoothing occurs only in horizontal direction
798    num_ax = len(scalar_grid.shape)
799    # Assume the last two axes represent the horizontal directions
800    sgma_seq = [sgma if i > num_ax - 3 else 0 for i in range(num_ax)]
801
802    # Compute smoothed field
803    return gaussian_filter(scalar_grid, sgma_seq, truncate=2 * np.sqrt(2))
804
805
806@exporter.export
807@preprocess_and_wrap(wrap_like='scalar_grid', match_unit=True, to_magnitude=True)
808def smooth_window(scalar_grid, window, passes=1, normalize_weights=True):
809    """Filter with an arbitrary window smoother.
810
811    Parameters
812    ----------
813    scalar_grid : array-like
814        N-dimensional scalar grid to be smoothed
815
816    window : ndarray
817        Window to use in smoothing. Can have dimension less than or equal to N. If
818        dimension less than N, the scalar grid will be smoothed along its trailing dimensions.
819        Shape along each dimension must be odd.
820
821    passes : int
822        The number of times to apply the filter to the grid. Defaults to 1.
823
824    normalize_weights : bool
825        If true, divide the values in window by the sum of all values in the window to obtain
826        the normalized smoothing weights. If false, use supplied values directly as the
827        weights.
828
829    Returns
830    -------
831    array-like
832        The filtered scalar grid
833
834    Notes
835    -----
836    This function can be applied multiple times to create a more smoothed field and will only
837    smooth the interior points, leaving the end points with their original values (this
838    function will leave an unsmoothed edge of size `(n - 1) / 2` for each `n` in the shape of
839    `window` around the data). If a masked value or NaN values exists in the array, it will
840    propagate to any point that uses that particular grid point in the smoothing calculation.
841    Applying the smoothing function multiple times will propagate NaNs further throughout the
842    domain.
843
844    See Also
845    --------
846    smooth_rectangular, smooth_circular, smooth_n_point, smooth_gaussian
847
848    """
849    def _pad(n):
850        # Return number of entries to pad given length along dimension.
851        return (n - 1) // 2
852
853    def _zero_to_none(x):
854        # Convert zero values to None, otherwise return what is given.
855        return x if x != 0 else None
856
857    def _offset(pad, k):
858        # Return padded slice offset by k entries
859        return slice(_zero_to_none(pad + k), _zero_to_none(-pad + k))
860
861    def _trailing_dims(indexer):
862        # Add ... to the front of an indexer, since we are working with trailing dimensions.
863        return (Ellipsis,) + tuple(indexer)
864
865    # Verify that shape in all dimensions is odd (need to have a neighborhood around a
866    # central point)
867    if any((size % 2 == 0) for size in window.shape):
868        raise ValueError('The shape of the smoothing window must be odd in all dimensions.')
869
870    # Optionally normalize the supplied weighting window
871    if normalize_weights:
872        weights = window / np.sum(window)
873    else:
874        weights = window
875
876    # Set indexes
877    # Inner index for the centered array elements that are affected by the smoothing
878    inner_full_index = _trailing_dims(_offset(_pad(n), 0) for n in weights.shape)
879    # Indexes to iterate over each weight
880    weight_indexes = tuple(product(*(range(n) for n in weights.shape)))
881
882    # Index for full array elements, offset by the weight index
883    def offset_full_index(weight_index):
884        return _trailing_dims(_offset(_pad(n), weight_index[i] - _pad(n))
885                              for i, n in enumerate(weights.shape))
886
887    # TODO: this is not lazy-loading/dask compatible, as it "densifies" the data
888    data = np.array(scalar_grid)
889    for _ in range(passes):
890        # Set values corresponding to smoothing weights by summing over each weight and
891        # applying offsets in needed dimensions
892        data[inner_full_index] = sum(weights[index] * data[offset_full_index(index)]
893                                     for index in weight_indexes)
894
895    return data
896
897
898@exporter.export
899def smooth_rectangular(scalar_grid, size, passes=1):
900    """Filter with a rectangular window smoother.
901
902    Parameters
903    ----------
904    scalar_grid : array-like
905        N-dimensional scalar grid to be smoothed
906
907    size : int or sequence of ints
908        Shape of rectangle along the trailing dimension(s) of the scalar grid
909
910    passes : int
911        The number of times to apply the filter to the grid. Defaults to 1.
912
913    Returns
914    -------
915    array-like
916        The filtered scalar grid
917
918    Notes
919    -----
920    This function can be applied multiple times to create a more smoothed field and will only
921    smooth the interior points, leaving the end points with their original values (this
922    function will leave an unsmoothed edge of size `(n - 1) / 2` for each `n` in `size` around
923    the data). If a masked value or NaN values exists in the array, it will propagate to any
924    point that uses that particular grid point in the smoothing calculation. Applying the
925    smoothing function multiple times will propagate NaNs further throughout the domain.
926
927    See Also
928    --------
929    smooth_window, smooth_circular, smooth_n_point, smooth_gaussian
930
931    """
932    return smooth_window(scalar_grid, np.ones(size), passes=passes)
933
934
935@exporter.export
936def smooth_circular(scalar_grid, radius, passes=1):
937    """Filter with a circular window smoother.
938
939    Parameters
940    ----------
941    scalar_grid : array-like
942        N-dimensional scalar grid to be smoothed. If more than two axes, smoothing is only
943        done along the last two.
944
945    radius : int
946        Radius of the circular smoothing window. The "diameter" of the circle (width of
947        smoothing window) is 2 * radius + 1 to provide a smoothing window with odd shape.
948
949    passes : int
950        The number of times to apply the filter to the grid. Defaults to 1.
951
952    Returns
953    -------
954    array-like
955        The filtered scalar grid
956
957    Notes
958    -----
959    This function can be applied multiple times to create a more smoothed field and will only
960    smooth the interior points, leaving the end points with their original values (this
961    function will leave an unsmoothed edge of size `radius` around the data). If a masked
962    value or NaN values exists in the array, it will propagate to any point that uses that
963    particular grid point in the smoothing calculation. Applying the smoothing function
964    multiple times will propagate NaNs further throughout the domain.
965
966    See Also
967    --------
968    smooth_window, smooth_rectangular, smooth_n_point, smooth_gaussian
969
970    """
971    # Generate the circle
972    size = 2 * radius + 1
973    x, y = np.mgrid[:size, :size]
974    distance = np.sqrt((x - radius) ** 2 + (y - radius) ** 2)
975    circle = distance <= radius
976
977    # Apply smoother
978    return smooth_window(scalar_grid, circle, passes=passes)
979
980
981@exporter.export
982def smooth_n_point(scalar_grid, n=5, passes=1):
983    """Filter with an n-point smoother.
984
985    Parameters
986    ----------
987    scalar_grid : array-like or `pint.Quantity`
988        N-dimensional scalar grid to be smoothed. If more than two axes, smoothing is only
989        done along the last two.
990
991    n: int
992        The number of points to use in smoothing, only valid inputs
993        are 5 and 9. Defaults to 5.
994
995    passes : int
996        The number of times to apply the filter to the grid. Defaults to 1.
997
998    Returns
999    -------
1000    array-like or `pint.Quantity`
1001        The filtered scalar grid
1002
1003    Notes
1004    -----
1005    This function is a close replication of the GEMPAK function SM5S and SM9S depending on the
1006    choice of the number of points to use for smoothing. This function can be applied multiple
1007    times to create a more smoothed field and will only smooth the interior points, leaving
1008    the end points with their original values (this function will leave an unsmoothed edge of
1009    size 1 around the data). If a masked value or NaN values exists in the array, it will
1010    propagate to any point that uses that particular grid point in the smoothing calculation.
1011    Applying the smoothing function multiple times will propagate NaNs further throughout the
1012    domain.
1013
1014    See Also
1015    --------
1016    smooth_window, smooth_rectangular, smooth_circular, smooth_gaussian
1017
1018    """
1019    if n == 9:
1020        weights = np.array([[0.0625, 0.125, 0.0625],
1021                            [0.125, 0.25, 0.125],
1022                            [0.0625, 0.125, 0.0625]])
1023    elif n == 5:
1024        weights = np.array([[0., 0.125, 0.],
1025                            [0.125, 0.5, 0.125],
1026                            [0., 0.125, 0.]])
1027    else:
1028        raise ValueError('The number of points to use in the smoothing '
1029                         'calculation must be either 5 or 9.')
1030
1031    return smooth_window(scalar_grid, window=weights, passes=passes, normalize_weights=False)
1032
1033
1034@exporter.export
1035@preprocess_and_wrap(wrap_like='altimeter_value')
1036@check_units('[pressure]', '[length]')
1037def altimeter_to_station_pressure(altimeter_value, height):
1038    r"""Convert the altimeter measurement to station pressure.
1039
1040    This function is useful for working with METARs since they do not provide
1041    altimeter values, but not sea-level pressure or station pressure.
1042    The following definitions of altimeter setting and station pressure
1043    are taken from [Smithsonian1951]_ Altimeter setting is the
1044    pressure value to which an aircraft altimeter scale is set so that it will
1045    indicate the altitude above mean sea-level of an aircraft on the ground at the
1046    location for which the value is determined. It assumes a standard atmosphere [NOAA1976]_.
1047    Station pressure is the atmospheric pressure at the designated station elevation.
1048    Finding the station pressure can be helpful for calculating sea-level pressure
1049    or other parameters.
1050
1051    Parameters
1052    ----------
1053    altimeter_value : `pint.Quantity`
1054        The altimeter setting value as defined by the METAR or other observation,
1055        which can be measured in either inches of mercury (in. Hg) or millibars (mb)
1056
1057    height: `pint.Quantity`
1058        Elevation of the station measuring pressure
1059
1060    Returns
1061    -------
1062    `pint.Quantity`
1063        The station pressure in hPa or in. Hg. Can be used to calculate sea-level
1064        pressure.
1065
1066    See Also
1067    --------
1068    altimeter_to_sea_level_pressure
1069
1070    Notes
1071    -----
1072    This function is implemented using the following equations from the
1073    Smithsonian Handbook (1951) p. 269
1074
1075    Equation 1:
1076     .. math:: A_{mb} = (p_{mb} - 0.3)F
1077
1078    Equation 3:
1079     .. math::  F = \left [1 + \left(\frac{p_{0}^n a}{T_{0}} \right)
1080                   \frac{H_{b}}{p_{1}^n} \right ] ^ \frac{1}{n}
1081
1082    Where,
1083
1084    :math:`p_{0}` = standard sea-level pressure = 1013.25 mb
1085
1086    :math:`p_{1} = p_{mb} - 0.3` when :math:`p_{0} = 1013.25 mb`
1087
1088    gamma = lapse rate in [NOAA1976]_ standard atmosphere below the isothermal layer
1089    :math:`6.5^{\circ}C. km.^{-1}`
1090
1091    :math:`t_{0}` = standard sea-level temperature 288 K
1092
1093    :math:`H_{b} =` station elevation in meters (elevation for which station pressure is given)
1094
1095    :math:`n = \frac{a R_{d}}{g} = 0.190284` where :math:`R_{d}` is the gas constant for dry
1096    air
1097
1098    And solving for :math:`p_{mb}` results in the equation below, which is used to
1099    calculate station pressure :math:`(p_{mb})`
1100
1101    .. math:: p_{mb} = \left [A_{mb} ^ n - \left (\frac{p_{0} a H_{b}}{T_0}
1102                       \right) \right] ^ \frac{1}{n} + 0.3
1103
1104    """
1105    # N-Value
1106    n = (mpconsts.Rd * gamma / mpconsts.g).to_base_units()
1107
1108    return ((altimeter_value ** n
1109             - ((p0.to(altimeter_value.units) ** n * gamma * height) / t0)) ** (1 / n)
1110            + units.Quantity(0.3, 'hPa'))
1111
1112
1113@exporter.export
1114@preprocess_and_wrap(wrap_like='altimeter_value')
1115@check_units('[pressure]', '[length]', '[temperature]')
1116def altimeter_to_sea_level_pressure(altimeter_value, height, temperature):
1117    r"""Convert the altimeter setting to sea-level pressure.
1118
1119    This function is useful for working with METARs since most provide
1120    altimeter values, but not sea-level pressure, which is often plotted
1121    on surface maps. The following definitions of altimeter setting, station pressure, and
1122    sea-level pressure are taken from [Smithsonian1951]_.
1123    Altimeter setting is the pressure value to which an aircraft altimeter scale
1124    is set so that it will indicate the altitude above mean sea-level of an aircraft
1125    on the ground at the location for which the value is determined. It assumes a standard
1126    atmosphere. Station pressure is the atmospheric pressure at the designated station
1127    elevation. Sea-level pressure is a pressure value obtained by the theoretical reduction
1128    of barometric pressure to sea level. It is assumed that atmosphere extends to sea level
1129    below the station and that the properties of the atmosphere are related to conditions
1130    observed at the station. This value is recorded by some surface observation stations,
1131    but not all. If the value is recorded, it can be found in the remarks section. Finding
1132    the sea-level pressure is helpful for plotting purposes and different calculations.
1133
1134    Parameters
1135    ----------
1136    altimeter_value : 'pint.Quantity'
1137        The altimeter setting value is defined by the METAR or other observation,
1138        with units of inches of mercury (in Hg) or millibars (hPa).
1139
1140    height  : 'pint.Quantity'
1141        Elevation of the station measuring pressure. Often times measured in meters
1142
1143    temperature : 'pint.Quantity'
1144        Temperature at the station
1145
1146    Returns
1147    -------
1148    'pint.Quantity'
1149        The sea-level pressure in hPa and makes pressure values easier to compare
1150        between different stations.
1151
1152    See Also
1153    --------
1154    altimeter_to_station_pressure
1155
1156    Notes
1157    -----
1158    This function is implemented using the following equations from Wallace and Hobbs (1977).
1159
1160    Equation 2.29:
1161     .. math::
1162       \Delta z = Z_{2} - Z_{1}
1163       = \frac{R_{d} \bar T_{v}}{g_0}ln\left(\frac{p_{1}}{p_{2}}\right)
1164       = \bar H ln \left (\frac {p_{1}}{p_{2}} \right)
1165
1166    Equation 2.31:
1167     .. math::
1168       p_{0} = p_{g}exp \left(\frac{Z_{g}}{\bar H} \right)
1169       = p_{g}exp \left(\frac{g_{0}Z_{g}}{R_{d}\bar T_{v}} \right)
1170
1171    Then by substituting :math:`\Delta_{Z}` for :math:`Z_{g}` in Equation 2.31:
1172     .. math:: p_{sealevel} = p_{station} exp\left(\frac{\Delta z}{H}\right)
1173
1174    where :math:`\Delta_{Z}` is the elevation in meters and :math:`H = \frac{R_{d}T}{g}`
1175
1176    """
1177    # Calculate the station pressure using function altimeter_to_station_pressure()
1178    psfc = altimeter_to_station_pressure(altimeter_value, height)
1179
1180    # Calculate the scale height
1181    h = mpconsts.Rd * temperature / mpconsts.g
1182
1183    return psfc * np.exp(height / h)
1184
1185
1186def _check_radians(value, max_radians=2 * np.pi):
1187    """Input validation of values that could be in degrees instead of radians.
1188
1189    Parameters
1190    ----------
1191    value : `pint.Quantity`
1192        Input value to check
1193
1194    max_radians : float
1195        Maximum absolute value of radians before warning
1196
1197    Returns
1198    -------
1199    `pint.Quantity`
1200        Input value
1201
1202    """
1203    with contextlib.suppress(AttributeError):
1204        value = value.to('radians').m
1205    if np.any(np.greater(np.abs(value), max_radians)):
1206        warnings.warn('Input over {} radians. '
1207                      'Ensure proper units are given.'.format(np.nanmax(max_radians)))
1208    return value
1209