add_pressure_to_height() - Code Metrics - Inspection of "Have get_upper_air_data bring in height data." - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#459)

unknown

created 2017-06-19 21:18 UTC

add_pressure_to_height() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
dl	0
loc	26
rs	8.8571
c	0
b	0
f	0

# Copyright (c) 2008-2015 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Contains a collection of basic calculations.

These include:

* wind components
* heat index
* windchill
"""

from __future__ import division

import numpy as np

from ..constants import g, omega, Rd
from ..package_tools import Exporter
from ..units import atleast_1d, check_units, masked_array, units

exporter = Exporter(globals())


@exporter.export
def get_wind_speed(u, v):
    r"""Compute the wind speed from u and v-components.

    Parameters
    ----------
    u : array_like
        Wind component in the X (East-West) direction
    v : array_like
        Wind component in the Y (North-South) direction

    Returns
    -------
    wind speed: array_like
        The speed of the wind

    See Also
    --------
    get_wind_components

    """
    speed = np.sqrt(u * u + v * v)
    return speed


@exporter.export
def get_wind_dir(u, v):
    r"""Compute the wind direction from u and v-components.

    Parameters
    ----------
    u : array_like
        Wind component in the X (East-West) direction
    v : array_like
        Wind component in the Y (North-South) direction

    Returns
    -------
    wind direction: `pint.Quantity`
        The direction of the wind, specified as the direction from
        which it is blowing, with 0 being North.

    See Also
    --------
    get_wind_components

    """
    wdir = 90. * units.deg - np.arctan2(-v, -u)
    origshape = wdir.shape
    wdir = atleast_1d(wdir)
    wdir[wdir < 0] += 360. * units.deg
    return wdir.reshape(origshape)


@exporter.export
def get_wind_components(speed, wdir):
    r"""Calculate the U, V wind vector components from the speed and direction.

    Parameters
    ----------
    speed : array_like
        The wind speed (magnitude)
    wdir : array_like
        The wind direction, specified as the direction from which the wind is
        blowing, with 0 being North.

    Returns
    -------
    u, v : tuple of array_like
        The wind components in the X (East-West) and Y (North-South)
        directions, respectively.

    See Also
    --------
    get_wind_speed
    get_wind_dir

    Examples
    --------
    >>> from metpy.units import units
    >>> metpy.calc.get_wind_components(10. * units('m/s'), 225. * units.deg)
    (<Quantity(7.071067811865475, 'meter / second')>,
     <Quantity(7.071067811865477, 'meter / second')>)

    """
    u = -speed * np.sin(wdir)
    v = -speed * np.cos(wdir)
    return u, v


@exporter.export
@check_units(temperature='[temperature]', speed='[speed]')
def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
    r"""Calculate the Wind Chill Temperature Index (WCTI).

    Calculates WCTI from the current temperature and wind speed using the formula
    outlined by the FCM [FCMR192003]_.

    Specifically, these formulas assume that wind speed is measured at
    10m.  If, instead, the speeds are measured at face level, the winds
    need to be multiplied by a factor of 1.5 (this can be done by specifying
    `face_level_winds` as `True`.)

    Parameters
    ----------
    temperature : `pint.Quantity`
        The air temperature
    speed : `pint.Quantity`
        The wind speed at 10m.  If instead the winds are at face level,
        `face_level_winds` should be set to `True` and the 1.5 multiplicative
        correction will be applied automatically.
    face_level_winds : bool, optional
        A flag indicating whether the wind speeds were measured at facial
        level instead of 10m, thus requiring a correction.  Defaults to
        `False`.
    mask_undefined : bool, optional
        A flag indicating whether a masked array should be returned with
        values where wind chill is undefined masked.  These are values where
        the temperature > 50F or wind speed <= 3 miles per hour. Defaults
        to `True`.

    Returns
    -------
    `pint.Quantity`
        The corresponding Wind Chill Temperature Index value(s)

    See Also
    --------
    heat_index

    """
    # Correct for lower height measurement of winds if necessary
    if face_level_winds:
        # No in-place so that we copy
        # noinspection PyAugmentAssignment
        speed = speed * 1.5

    temp_limit, speed_limit = 10. * units.degC, 3 * units.mph
    speed_factor = speed.to('km/hr').magnitude ** 0.16
    wcti = units.Quantity((0.6215 + 0.3965 * speed_factor) * temperature.to('degC').magnitude -
                          11.37 * speed_factor + 13.12, units.degC).to(temperature.units)

    # See if we need to mask any undefined values
    if mask_undefined:
        mask = np.array((temperature > temp_limit) | (speed <= speed_limit))
        if mask.any():
            wcti = masked_array(wcti, mask=mask)

    return wcti


@exporter.export
@check_units('[temperature]')
def heat_index(temperature, rh, mask_undefined=True):
    r"""Calculate the Heat Index from the current temperature and relative humidity.

    The implementation uses the formula outlined in [Rothfusz1990]_. This equation is a
    multi-variable least-squares regression of the values obtained in [Steadman1979]_.

    Parameters
    ----------
    temperature : `pint.Quantity`
        Air temperature
    rh : array_like
        The relative humidity expressed as a unitless ratio in the range [0, 1].
        Can also pass a percentage if proper units are attached.

    Returns
    -------
    `pint.Quantity`
        The corresponding Heat Index value(s)

    Other Parameters
    ----------------
    mask_undefined : bool, optional
        A flag indicating whether a masked array should be returned with
        values where heat index is undefined masked.  These are values where
        the temperature < 80F or relative humidity < 40 percent. Defaults
        to `True`.

    See Also
    --------
    windchill

    """
    delta = temperature - 0. * units.degF
    rh2 = rh * rh
    delta2 = delta * delta

    # Calculate the Heat Index -- constants converted for RH in [0, 1]
    hi = (-42.379 * units.degF + 2.04901523 * delta +
          1014.333127 * units.delta_degF * rh - 22.475541 * delta * rh -
          6.83783e-3 / units.delta_degF * delta2 - 5.481717e2 * units.delta_degF * rh2 +
          1.22874e-1 / units.delta_degF * delta2 * rh + 8.5282 * delta * rh2 -
          1.99e-2 / units.delta_degF * delta2 * rh2)

    # See if we need to mask any undefined values
    if mask_undefined:
        mask = np.array((temperature < 80. * units.degF) | (rh < 40 * units.percent))
        if mask.any():
            hi = masked_array(hi, mask=mask)

    return hi


@exporter.export
@check_units('[pressure]')
def pressure_to_height_std(pressure):
    r"""Convert pressure data to heights using the U.S. standard atmosphere.

    The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61.

    Parameters
    ----------
    pressure : `pint.Quantity`
        Atmospheric pressure

    Returns
    -------
    `pint.Quantity`
        The corresponding height value(s)

    Notes
    -----
    .. math:: Z = \frac{T_0}{\Gamma}[1-\frac{p}{p_0}^\frac{R\Gamma}{g}]

    """
    t0 = 288. * units.kelvin
    gamma = 6.5 * units('K/km')
    p0 = 1013.25 * units.mbar
    return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(Rd * gamma / g))


@exporter.export
@check_units('[length]')
def height_to_pressure_std(height):
    r"""Convert height data to pressures using the U.S. standard atmosphere.

    The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61.

    Parameters
    ----------
    height : `pint.Quantity`
        Atmospheric height

    Returns
    -------
    `pint.Quantity`
        The corresponding pressure value(s)

    Notes
    -----
    .. math:: p = p_0 e^{\frac{g}{R \Gamma} \text{ln}(1-\frac{Z \Gamma}{T_0})}

    """
    t0 = 288. * units.kelvin
    gamma = 6.5 * units('K/km')
    p0 = 1013.25 * units.mbar
    return p0 * (1 - (gamma / t0) * height) ** (g / (Rd * gamma))


@exporter.export
def coriolis_parameter(latitude):
    r"""Calculate the coriolis parameter at each point.

    The implementation uses the formula outlined in [Hobbs1977]_ pg.370-371.

    Parameters
    ----------
    latitude : array_like
        Latitude at each point

    Returns
    -------
    `pint.Quantity`
        The corresponding coriolis force at each point

    """
    return 2. * omega * np.sin(latitude)


@exporter.export
@check_units('[pressure]', '[length]')
def add_height_to_pressure(pressure, height):
    r"""Calculate the pressure at a certain height above another pressure level.

    This assumes a standard atmosphere.

    Parameters
    ----------
    pressure : `pint.Quantity`
        Pressure level
    height : `pint.Quantity`
        Height above a pressure level

    Returns
    -------
    `pint.Quantity`
        The corresponding pressure value for the height above the pressure level

    See Also
    -----
    pressure_to_height_std, height_to_pressure_std, add_pressure_to_height

    """
    pressure_level_height = pressure_to_height_std(pressure)
    return height_to_pressure_std(pressure_level_height + height)


@exporter.export
@check_units('[length]', '[pressure]')
def add_pressure_to_height(height, pressure):
    r"""Calculate the height at a certain pressure above another height.

    This assumes a standard atmosphere.

    Parameters
    ----------
    height : `pint.Quantity`
        Height level
    pressure : `pint.Quantity`
        Pressure above height level

    Returns
    -------
    `pint.Quantity`
        The corresponding height value for the pressure above the height level

    See Also
    -----
    pressure_to_height_std, height_to_pressure_std, add_height_to_pressure

    """
    pressure_at_height = height_to_pressure_std(height)
    return pressure_to_height_std(pressure_at_height - pressure)


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

Unidata / MetPy

Pull Request — master (#459)

add_pressure_to_height() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like