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