el() - Code Metrics - Inspection of "CAPE and CIN calculations" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#322)

unknown

created 2017-02-15 17:31 UTC

el() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

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

# 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 thermodynamic calculations."""

from __future__ import division

import numpy as np
import scipy.integrate as si

from .tools import find_intersections
from ..constants import Cp_d, epsilon, kappa, Lv, P0, Rd
from ..package_tools import Exporter
from ..units import atleast_1d, concatenate, units

exporter = Exporter(globals())

sat_pressure_0c = 6.112 * units.millibar


@exporter.export
def potential_temperature(pressure, temperature):
    r"""Calculate the potential temperature.

    Uses the Poisson equation to calculation the potential temperature
    given `pressure` and `temperature`.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The total atmospheric pressure
    temperature : `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The potential temperature corresponding to the the temperature and
        pressure.

    See Also
    --------
    dry_lapse

    Notes
    -----
    Formula:

    .. math:: \Theta = T (P_0 / P)^\kappa

    Examples
    --------
    >>> from metpy.units import units
    >>> metpy.calc.potential_temperature(800. * units.mbar, 273. * units.kelvin)
    290.9814150577374
    """
    return temperature * (P0 / pressure).to('dimensionless')**kappa


@exporter.export
def dry_lapse(pressure, temperature):
    r"""Calculate the temperature at a level assuming only dry processes.

    This function lifts a parcel starting at `temperature`, conserving
    potential temperature. The starting pressure should be the first item in
    the `pressure` array.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest
    temperature : `pint.Quantity`
        The starting temperature

    Returns
    -------
    `pint.Quantity`
       The resulting parcel temperature at levels given by `pressure`

    See Also
    --------
    moist_lapse : Calculate parcel temperature assuming liquid saturation
                  processes
    parcel_profile : Calculate complete parcel profile
    potential_temperature
    """
    return temperature * (pressure / pressure[0])**kappa


@exporter.export
def moist_lapse(pressure, temperature):
    r"""Calculate the temperature at a level assuming liquid saturation processes.

    This function lifts a parcel starting at `temperature`. The starting
    pressure should be the first item in the `pressure` array. Essentially,
    this function is calculating moist pseudo-adiabats.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest
    temperature : `pint.Quantity`
        The starting temperature

    Returns
    -------
    `pint.Quantity`
       The temperature corresponding to the the starting temperature and
       pressure levels.

    See Also
    --------
    dry_lapse : Calculate parcel temperature assuming dry adiabatic processes
    parcel_profile : Calculate complete parcel profile

    Notes
    -----
    This function is implemented by integrating the following differential
    equation:

    .. math:: \frac{dT}{dP} = \frac{1}{P} \frac{R_d T + L_v r_s}
                                {C_{pd} + \frac{L_v^2 r_s \epsilon}{R_d T^2}}

    This equation comes from [1]_.

    References
    ----------
    .. [1] Bakhshaii, A. and R. Stull, 2013: Saturated Pseudoadiabats--A
           Noniterative Approximation. J. Appl. Meteor. Clim., 52, 5-15.
    """
    def dt(t, p):
        t = units.Quantity(t, temperature.units)
        p = units.Quantity(p, pressure.units)
        rs = saturation_mixing_ratio(p, t)
        frac = ((Rd * t + Lv * rs) /
                (Cp_d + (Lv * Lv * rs * epsilon / (Rd * t * t)))).to('kelvin')
        return frac / p
    return units.Quantity(si.odeint(dt, atleast_1d(temperature).squeeze(),
                                    pressure.squeeze()).T.squeeze(), temperature.units)


@exporter.export
def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-2):
    r"""Calculate the lifted condensation level (LCL) using from the starting point.

    The starting state for the parcel is defined by `temperature`, `dewpt`,
    and `pressure`.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The starting atmospheric pressure
    temperature : `pint.Quantity`
        The starting temperature
    dewpt : `pint.Quantity`
        The starting dew point

    Returns
    -------
    `(pint.Quantity, pint.Quantity)`
        The LCL pressure and temperature

    Other Parameters
    ----------------
    max_iters : int, optional
        The maximum number of iterations to use in calculation, defaults to 50.
    eps : float, optional
        The desired absolute error in the calculated value, defaults to 1e-2.

    See Also
    --------
    parcel_profile

    Notes
    -----
    This function is implemented using an iterative approach to solve for the
    LCL. The basic algorithm is:

    1. Find the dew point from the LCL pressure and starting mixing ratio
    2. Find the LCL pressure from the starting temperature and dewpoint
    3. Iterate until convergence

    The function is guaranteed to finish by virtue of the `max_iters` counter.
    """
    w = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
    new_p = p = pressure
    eps = units.Quantity(eps, p.units)
    while max_iters:
        td = dewpoint(vapor_pressure(p, w))
        new_p = pressure * (td / temperature) ** (1. / kappa)
        if np.abs(new_p - p).max() < eps:
            break
        p = new_p
        max_iters -= 1

    return (new_p, td)


@exporter.export
def lfc(pressure, temperature, dewpt):
    r"""Calculate the level of free convection (LFC).

    This works by finding the first intersection of the ideal parcel path and
    the measured parcel temperature.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure
    temperature : `pint.Quantity`
        The temperature at the levels given by `pressure`
    dewpt : `pint.Quantity`
        The dew point at the levels given by `pressure`

    Returns
    -------
    `pint.Quantity`
        The LFC

    See Also
    --------
    parcel_profile
    """
    ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
    x, y = find_intersections(pressure, ideal_profile, temperature)
    return x[0], y[0]


@exporter.export
def el(pressure, temperature, dewpt):
    r"""Calculate the equilibrium level (EL).

    This works by finding the last intersection of the ideal parcel path and
    the measured parcel temperature. If there is one or fewer intersections, there is
    no equilibrium level.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure
    temperature : `pint.Quantity`
        The temperature at the levels given by `pressure`
    dewpt : `pint.Quantity`
        The dew point at the levels given by `pressure`

    Returns
    -------
    `pint.Quantity, pint.Quantity`
        The EL pressure and temperature

    See Also
    --------
    parcel_profile
    """
    ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
    x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:])

    # If there is only one intersection, it's the LFC and we return None.
    if len(x) <= 1:
        return None, None

    else:
        return x[-1], y[-1]


@exporter.export
def parcel_profile(pressure, temperature, dewpt):
    r"""Calculate the profile a parcel takes through the atmosphere.

    The parcel starts at `temperature`, and `dewpt`, lifted up
    dry adiabatically to the LCL, and then moist adiabatically from there.
    `pressure` specifies the pressure levels for the profile.

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest. The first entry should be the starting
        point pressure.
    temperature : `pint.Quantity`
        The starting temperature
    dewpt : `pint.Quantity`
        The starting dew point

    Returns
    -------
    `pint.Quantity`
        The parcel temperatures at the specified pressure levels.

    See Also
    --------
    lcl, moist_lapse, dry_lapse
    """
    # Find the LCL
    l = lcl(pressure[0], temperature, dewpt)[0].to(pressure.units)

    # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the
    # LCL is included in the levels. It's slightly redundant in that case, but simplifies
    # the logic for removing it later.
    press_lower = concatenate((pressure[pressure >= l], l))
    t1 = dry_lapse(press_lower, temperature)

    # Find moist pseudo-adiabatic profile starting at the LCL
    press_upper = concatenate((l, pressure[pressure < l]))
    t2 = moist_lapse(press_upper, t1[-1]).to(t1.units)

    # Return LCL *without* the LCL point
    return concatenate((t1[:-1], t2[1:]))


@exporter.export
def virtual_temperature(pressure, temperature, dewpt):
    r"""Calculate virtual temperature.

    Description

    Parameters
    ----------
    pressure : `pint.Quantity`
        The atmospheric pressure level(s) of interest. The first entry should be the starting
        point pressure.
    temperature : `pint.Quantity`
        The starting temperature
    dewpt : `pint.Quantity`
        The starting dew point

    Returns
    -------
    `pint.Quantity`
        The virtual temperature.

    See Also
    --------
    mixing_ratio, saturation_vapor_pressure
    """
    Rv = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
    return temperature * (1 + Rv / epsilon) / (1 + Rv)


@exporter.export
def cape(pressure, temperature, dewpt, temperature_parcel, virtual_temp=False):

    if virtual_temp:
        temperature = virtual_temperature(pressure, temperature, dewpt)
        temperature_parcel = virtual_temperature(pressure, temperature_parcel, dewpt)

    x = np.copy(pressure)
    y = temperature_parcel - temperature

    # Find and append crossings to the data
    crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
    x = np.concatenate((x, crossings[0]))
    y = np.concatenate((y, crossings[1]))

    # Resort so that data is in pressure order
    sort_idx = np.argsort(x)
    x = x[sort_idx]
    y = y[sort_idx]

    # Remove duplicate data points if there are any
    keep_idx = np.ediff1d(x, to_end=1) >= 1
    x = x[keep_idx]
    y = y[keep_idx]

    # Clip out negative area (temperature parcel < temperature environment)
    y[y <= 0] = 0

    # Only use data between the LFC and EL for calculation
    lfc_pressure = lfc(pressure, temperature, dewpt)[0].magnitude
    el_pressure = el(pressure, temperature, dewpt)[0].magnitude
    p_mask = (x <= lfc_pressure) & (x >= el_pressure)
    x = x[p_mask]
    y = y[p_mask]

    return Rd * np.trapz(y, np.log(x))


@exporter.export
def vapor_pressure(pressure, mixing):
    r"""Calculate water vapor (partial) pressure.

    Given total `pressure` and water vapor `mixing` ratio, calculates the
    partial pressure of water vapor.

    Parameters
    ----------
    pressure : `pint.Quantity`
        total atmospheric pressure
    mixing : `pint.Quantity`
        dimensionless mass mixing ratio

    Returns
    -------
    `pint.Quantity`
        The ambient water vapor (partial) pressure in the same units as
        `pressure`.

    Notes
    -----
    This function is a straightforward implementation of the equation given in many places,
    such as [2]_:

    .. math:: e = p \frac{r}{r + \epsilon}

    References
    ----------
    .. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 71.

    See Also
    --------
    saturation_vapor_pressure, dewpoint
    """
    return pressure * mixing / (epsilon + mixing)


@exporter.export
def saturation_vapor_pressure(temperature):
    r"""Calculate the saturation water vapor (partial) pressure.

    Parameters
    ----------
    temperature : `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The saturation water vapor (partial) pressure

    See Also
    --------
    vapor_pressure, dewpoint

    Notes
    -----
    Instead of temperature, dewpoint may be used in order to calculate
    the actual (ambient) water vapor (partial) pressure.

    The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:

    .. math:: 6.112 e^\frac{17.67T}{T + 243.5}

    References
    ----------
    .. [3] Bolton, D., 1980: The Computation of Equivalent Potential
           Temperature. Mon. Wea. Rev., 108, 1046-1053.
    """
    # Converted from original in terms of C to use kelvin. Using raw absolute values of C in
    # a formula plays havoc with units support.
    return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
                                    (temperature - 29.65 * units.kelvin))


@exporter.export
def dewpoint_rh(temperature, rh):
    r"""Calculate the ambient dewpoint given air temperature and relative humidity.

    Parameters
    ----------
    temperature : `pint.Quantity`
        Air temperature
    rh : `pint.Quantity`
        Relative humidity expressed as a ratio in the range [0, 1]

    Returns
    -------
    `pint.Quantity`
        The dew point temperature

    See Also
    --------
    dewpoint, saturation_vapor_pressure
    """
    return dewpoint(rh * saturation_vapor_pressure(temperature))


@exporter.export
def dewpoint(e):
    r"""Calculate the ambient dewpoint given the vapor pressure.

    Parameters
    ----------
    e : `pint.Quantity`
        Water vapor partial pressure

    Returns
    -------
    `pint.Quantity`
        Dew point temperature

    See Also
    --------
    dewpoint_rh, saturation_vapor_pressure, vapor_pressure

    Notes
    -----
    This function inverts the Bolton 1980 [4]_ formula for saturation vapor
    pressure to instead calculate the temperature. This yield the following
    formula for dewpoint in degrees Celsius:

    .. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}

    References
    ----------
    .. [4] Bolton, D., 1980: The Computation of Equivalent Potential
           Temperature. Mon. Wea. Rev., 108, 1046-1053.
    """
    val = np.log(e / sat_pressure_0c)
    return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)


@exporter.export
def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
    r"""Calculate the mixing ratio of a gas.

    This calculates mixing ratio given its partial pressure and the total pressure of
    the air. There are no required units for the input arrays, other than that
    they have the same units.

    Parameters
    ----------
    part_press : `pint.Quantity`
        Partial pressure of the constituent gas
    tot_press : `pint.Quantity`
        Total air pressure
    molecular_weight_ratio : `pint.Quantity` or float, optional
        The ratio of the molecular weight of the constituent gas to that assumed
        for air. Defaults to the ratio for water vapor to dry air
        (:math:`\epsilon\approx0.622`).

    Returns
    -------
    `pint.Quantity`
        The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)

    Notes
    -----
    This function is a straightforward implementation of the equation given in many places,
    such as [5]_:

    .. math:: r = \epsilon \frac{e}{p - e}

    References
    ----------
    .. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 73.

    See Also
    --------
    saturation_mixing_ratio, vapor_pressure
    """
    return (molecular_weight_ratio * part_press / (tot_press - part_press)).to("dimensionless")


@exporter.export
def saturation_mixing_ratio(tot_press, temperature):
    r"""Calculate the saturation mixing ratio of water vapor.

    This calculation is given total pressure and the temperature. The implementation
    uses the formula outlined in [6]_.

    Parameters
    ----------
    tot_press: `pint.Quantity`
        Total atmospheric pressure
    temperature: `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The saturation mixing ratio, dimensionless

    References
    ----------
    .. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 73.
    """
    return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)


@exporter.export
def equivalent_potential_temperature(pressure, temperature):
    r"""Calculate equivalent potential temperature.

    This calculation must be given an air parcel's pressure and temperature.
    The implementation uses the formula outlined in [7]_.

    Parameters
    ----------
    pressure: `pint.Quantity`
        Total atmospheric pressure
    temperature: `pint.Quantity`
        The temperature

    Returns
    -------
    `pint.Quantity`
        The corresponding equivalent potential temperature of the parcel

    Notes
    -----
    .. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}

    References
    ----------
    .. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
           Survey. 78-79.
    """
    pottemp = potential_temperature(pressure, temperature)
    smixr = saturation_mixing_ratio(pressure, temperature)
    return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))


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 thermodynamic calculations."""
5
6			from __future__ import division
7
8			import numpy as np
9			import scipy.integrate as si
10
11			from .tools import find_intersections
12			from ..constants import Cp_d, epsilon, kappa, Lv, P0, Rd
13			from ..package_tools import Exporter
14			from ..units import atleast_1d, concatenate, units
15
16			exporter = Exporter(globals())
17
18			sat_pressure_0c = 6.112 * units.millibar
19
20
21			@exporter.export
22			def potential_temperature(pressure, temperature):
23			r"""Calculate the potential temperature.
24
25			Uses the Poisson equation to calculation the potential temperature
26			given `pressure` and `temperature`.
27
28			Parameters
29			----------
30			pressure : `pint.Quantity`
31			The total atmospheric pressure
32			temperature : `pint.Quantity`
33			The temperature
34
35			Returns
36			-------
37			`pint.Quantity`
38			The potential temperature corresponding to the the temperature and
39			pressure.
40
41			See Also
42			--------
43			dry_lapse
44
45			Notes
46			-----
47			Formula:
48
49			.. math:: \Theta = T (P_0 / P)^\kappa
50
51			Examples
52			--------
53			>>> from metpy.units import units
54			>>> metpy.calc.potential_temperature(800. * units.mbar, 273. * units.kelvin)
55			290.9814150577374
56			"""
57			return temperature * (P0 / pressure).to('dimensionless')**kappa
58
59
60			@exporter.export
61			def dry_lapse(pressure, temperature):
62			r"""Calculate the temperature at a level assuming only dry processes.
63
64			This function lifts a parcel starting at `temperature`, conserving
65			potential temperature. The starting pressure should be the first item in
66			the `pressure` array.
67
68			Parameters
69			----------
70			pressure : `pint.Quantity`
71			The atmospheric pressure level(s) of interest
72			temperature : `pint.Quantity`
73			The starting temperature
74
75			Returns
76			-------
77			`pint.Quantity`
78			The resulting parcel temperature at levels given by `pressure`
79
80			See Also
81			--------
82			moist_lapse : Calculate parcel temperature assuming liquid saturation
83			processes
84			parcel_profile : Calculate complete parcel profile
85			potential_temperature
86			"""
87			return temperature * (pressure / pressure[0])**kappa
88
89
90			@exporter.export
91			def moist_lapse(pressure, temperature):
92			r"""Calculate the temperature at a level assuming liquid saturation processes.
93
94			This function lifts a parcel starting at `temperature`. The starting
95			pressure should be the first item in the `pressure` array. Essentially,
96			this function is calculating moist pseudo-adiabats.
97
98			Parameters
99			----------
100			pressure : `pint.Quantity`
101			The atmospheric pressure level(s) of interest
102			temperature : `pint.Quantity`
103			The starting temperature
104
105			Returns
106			-------
107			`pint.Quantity`
108			The temperature corresponding to the the starting temperature and
109			pressure levels.
110
111			See Also
112			--------
113			dry_lapse : Calculate parcel temperature assuming dry adiabatic processes
114			parcel_profile : Calculate complete parcel profile
115
116			Notes
117			-----
118			This function is implemented by integrating the following differential
119			equation:
120
121			.. math:: \frac{dT}{dP} = \frac{1}{P} \frac{R_d T + L_v r_s}
122			{C_{pd} + \frac{L_v^2 r_s \epsilon}{R_d T^2}}
123
124			This equation comes from [1]_.
125
126			References
127			----------
128			.. [1] Bakhshaii, A. and R. Stull, 2013: Saturated Pseudoadiabats--A
129			Noniterative Approximation. J. Appl. Meteor. Clim., 52, 5-15.
130			"""
131			def dt(t, p):
132			t = units.Quantity(t, temperature.units)
133			p = units.Quantity(p, pressure.units)
134			rs = saturation_mixing_ratio(p, t)
135			frac = ((Rd * t + Lv * rs) /
136			(Cp_d + (Lv * Lv * rs * epsilon / (Rd * t * t)))).to('kelvin')
137			return frac / p
138			return units.Quantity(si.odeint(dt, atleast_1d(temperature).squeeze(),
139			pressure.squeeze()).T.squeeze(), temperature.units)
140
141
142			@exporter.export
143			def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-2):
144			r"""Calculate the lifted condensation level (LCL) using from the starting point.
145
146			The starting state for the parcel is defined by `temperature`, `dewpt`,
147			and `pressure`.
148
149			Parameters
150			----------
151			pressure : `pint.Quantity`
152			The starting atmospheric pressure
153			temperature : `pint.Quantity`
154			The starting temperature
155			dewpt : `pint.Quantity`
156			The starting dew point
157
158			Returns
159			-------
160			`(pint.Quantity, pint.Quantity)`
161			The LCL pressure and temperature
162
163			Other Parameters
164			----------------
165			max_iters : int, optional
166			The maximum number of iterations to use in calculation, defaults to 50.
167			eps : float, optional
168			The desired absolute error in the calculated value, defaults to 1e-2.
169
170			See Also
171			--------
172			parcel_profile
173
174			Notes
175			-----
176			This function is implemented using an iterative approach to solve for the
177			LCL. The basic algorithm is:
178
179			1. Find the dew point from the LCL pressure and starting mixing ratio
180			2. Find the LCL pressure from the starting temperature and dewpoint
181			3. Iterate until convergence
182
183			The function is guaranteed to finish by virtue of the `max_iters` counter.
184			"""
185			w = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
186			new_p = p = pressure
187			eps = units.Quantity(eps, p.units)
188			while max_iters:
189			td = dewpoint(vapor_pressure(p, w))
190			new_p = pressure * (td / temperature) ** (1. / kappa)
191			if np.abs(new_p - p).max() < eps:
192			break
193			p = new_p
194			max_iters -= 1
195
196			return (new_p, td)
197
198
199			@exporter.export
200			def lfc(pressure, temperature, dewpt):
201			r"""Calculate the level of free convection (LFC).
202
203			This works by finding the first intersection of the ideal parcel path and
204			the measured parcel temperature.
205
206			Parameters
207			----------
208			pressure : `pint.Quantity`
209			The atmospheric pressure
210			temperature : `pint.Quantity`
211			The temperature at the levels given by `pressure`
212			dewpt : `pint.Quantity`
213			The dew point at the levels given by `pressure`
214
215			Returns
216			-------
217			`pint.Quantity`
218			The LFC
219
220			See Also
221			--------
222			parcel_profile
223			"""
224			ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
225			x, y = find_intersections(pressure, ideal_profile, temperature)
226			return x[0], y[0]
227
228
229			@exporter.export
230			def el(pressure, temperature, dewpt):
231			r"""Calculate the equilibrium level (EL).
232
233			This works by finding the last intersection of the ideal parcel path and
234			the measured parcel temperature. If there is one or fewer intersections, there is
235			no equilibrium level.
236
237			Parameters
238			----------
239			pressure : `pint.Quantity`
240			The atmospheric pressure
241			temperature : `pint.Quantity`
242			The temperature at the levels given by `pressure`
243			dewpt : `pint.Quantity`
244			The dew point at the levels given by `pressure`
245
246			Returns
247			-------
248			`pint.Quantity, pint.Quantity`
249			The EL pressure and temperature
250
251			See Also
252			--------
253			parcel_profile
254			"""
255			ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC')
256			x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:])
257
258			# If there is only one intersection, it's the LFC and we return None.
259			if len(x) <= 1:
260			return None, None
261
262			else:
263			return x[-1], y[-1]
264
265
266			@exporter.export
267			def parcel_profile(pressure, temperature, dewpt):
268			r"""Calculate the profile a parcel takes through the atmosphere.
269
270			The parcel starts at `temperature`, and `dewpt`, lifted up
271			dry adiabatically to the LCL, and then moist adiabatically from there.
272			`pressure` specifies the pressure levels for the profile.
273
274			Parameters
275			----------
276			pressure : `pint.Quantity`
277			The atmospheric pressure level(s) of interest. The first entry should be the starting
278			point pressure.
279			temperature : `pint.Quantity`
280			The starting temperature
281			dewpt : `pint.Quantity`
282			The starting dew point
283
284			Returns
285			-------
286			`pint.Quantity`
287			The parcel temperatures at the specified pressure levels.
288
289			See Also
290			--------
291			lcl, moist_lapse, dry_lapse
292			"""
293			# Find the LCL
294			l = lcl(pressure[0], temperature, dewpt)[0].to(pressure.units)
295
296			# Find the dry adiabatic profile, including the LCL. We need >= the LCL in case the
297			# LCL is included in the levels. It's slightly redundant in that case, but simplifies
298			# the logic for removing it later.
299			press_lower = concatenate((pressure[pressure >= l], l))
300			t1 = dry_lapse(press_lower, temperature)
301
302			# Find moist pseudo-adiabatic profile starting at the LCL
303			press_upper = concatenate((l, pressure[pressure < l]))
304			t2 = moist_lapse(press_upper, t1[-1]).to(t1.units)
305
306			# Return LCL without the LCL point
307			return concatenate((t1[:-1], t2[1:]))
308
309
310			@exporter.export
311			def virtual_temperature(pressure, temperature, dewpt):
312			r"""Calculate virtual temperature.
313
314			Description
315
316			Parameters
317			----------
318			pressure : `pint.Quantity`
319			The atmospheric pressure level(s) of interest. The first entry should be the starting
320			point pressure.
321			temperature : `pint.Quantity`
322			The starting temperature
323			dewpt : `pint.Quantity`
324			The starting dew point
325
326			Returns
327			-------
328			`pint.Quantity`
329			The virtual temperature.
330
331			See Also
332			--------
333			mixing_ratio, saturation_vapor_pressure
334			"""
335			Rv = mixing_ratio(saturation_vapor_pressure(dewpt), pressure)
336			return temperature * (1 + Rv / epsilon) / (1 + Rv)
337
338
339			@exporter.export
340			def cape(pressure, temperature, dewpt, temperature_parcel, virtual_temp=False):
341
342			if virtual_temp:
343			temperature = virtual_temperature(pressure, temperature, dewpt)
344			temperature_parcel = virtual_temperature(pressure, temperature_parcel, dewpt)
345
346			x = np.copy(pressure)
347			y = temperature_parcel - temperature
348
349			# Find and append crossings to the data
350			crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
351			x = np.concatenate((x, crossings[0]))
352			y = np.concatenate((y, crossings[1]))
353
354			# Resort so that data is in pressure order
355			sort_idx = np.argsort(x)
356			x = x[sort_idx]
357			y = y[sort_idx]
358
359			# Remove duplicate data points if there are any
360			keep_idx = np.ediff1d(x, to_end=1) >= 1
361			x = x[keep_idx]
362			y = y[keep_idx]
363
364			# Clip out negative area (temperature parcel < temperature environment)
365			y[y <= 0] = 0
366
367			# Only use data between the LFC and EL for calculation
368			lfc_pressure = lfc(pressure, temperature, dewpt)[0].magnitude
369			el_pressure = el(pressure, temperature, dewpt)[0].magnitude
370			p_mask = (x <= lfc_pressure) & (x >= el_pressure)
371			x = x[p_mask]
372			y = y[p_mask]
373
374			return Rd * np.trapz(y, np.log(x))
375
376
377			@exporter.export
378			def vapor_pressure(pressure, mixing):
379			r"""Calculate water vapor (partial) pressure.
380
381			Given total `pressure` and water vapor `mixing` ratio, calculates the
382			partial pressure of water vapor.
383
384			Parameters
385			----------
386			pressure : `pint.Quantity`
387			total atmospheric pressure
388			mixing : `pint.Quantity`
389			dimensionless mass mixing ratio
390
391			Returns
392			-------
393			`pint.Quantity`
394			The ambient water vapor (partial) pressure in the same units as
395			`pressure`.
396
397			Notes
398			-----
399			This function is a straightforward implementation of the equation given in many places,
400			such as [2]_:
401
402			.. math:: e = p \frac{r}{r + \epsilon}
403
404			References
405			----------
406			.. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
407			Survey. 71.
408
409			See Also
410			--------
411			saturation_vapor_pressure, dewpoint
412			"""
413			return pressure * mixing / (epsilon + mixing)
414
415
416			@exporter.export
417			def saturation_vapor_pressure(temperature):
418			r"""Calculate the saturation water vapor (partial) pressure.
419
420			Parameters
421			----------
422			temperature : `pint.Quantity`
423			The temperature
424
425			Returns
426			-------
427			`pint.Quantity`
428			The saturation water vapor (partial) pressure
429
430			See Also
431			--------
432			vapor_pressure, dewpoint
433
434			Notes
435			-----
436			Instead of temperature, dewpoint may be used in order to calculate
437			the actual (ambient) water vapor (partial) pressure.
438
439			The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:
440
441			.. math:: 6.112 e^\frac{17.67T}{T + 243.5}
442
443			References
444			----------
445			.. [3] Bolton, D., 1980: The Computation of Equivalent Potential
446			Temperature. Mon. Wea. Rev., 108, 1046-1053.
447			"""
448			# Converted from original in terms of C to use kelvin. Using raw absolute values of C in
449			# a formula plays havoc with units support.
450			return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
451			(temperature - 29.65 * units.kelvin))
452
453
454			@exporter.export
455			def dewpoint_rh(temperature, rh):
456			r"""Calculate the ambient dewpoint given air temperature and relative humidity.
457
458			Parameters
459			----------
460			temperature : `pint.Quantity`
461			Air temperature
462			rh : `pint.Quantity`
463			Relative humidity expressed as a ratio in the range [0, 1]
464
465			Returns
466			-------
467			`pint.Quantity`
468			The dew point temperature
469
470			See Also
471			--------
472			dewpoint, saturation_vapor_pressure
473			"""
474			return dewpoint(rh * saturation_vapor_pressure(temperature))
475
476
477			@exporter.export
478			def dewpoint(e):
479			r"""Calculate the ambient dewpoint given the vapor pressure.
480
481			Parameters
482			----------
483			e : `pint.Quantity`
484			Water vapor partial pressure
485
486			Returns
487			-------
488			`pint.Quantity`
489			Dew point temperature
490
491			See Also
492			--------
493			dewpoint_rh, saturation_vapor_pressure, vapor_pressure
494
495			Notes
496			-----
497			This function inverts the Bolton 1980 [4]_ formula for saturation vapor
498			pressure to instead calculate the temperature. This yield the following
499			formula for dewpoint in degrees Celsius:
500
501			.. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}
502
503			References
504			----------
505			.. [4] Bolton, D., 1980: The Computation of Equivalent Potential
506			Temperature. Mon. Wea. Rev., 108, 1046-1053.
507			"""
508			val = np.log(e / sat_pressure_0c)
509			return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
510
511
512			@exporter.export
513			def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
514			r"""Calculate the mixing ratio of a gas.
515
516			This calculates mixing ratio given its partial pressure and the total pressure of
517			the air. There are no required units for the input arrays, other than that
518			they have the same units.
519
520			Parameters
521			----------
522			part_press : `pint.Quantity`
523			Partial pressure of the constituent gas
524			tot_press : `pint.Quantity`
525			Total air pressure
526			molecular_weight_ratio : `pint.Quantity` or float, optional
527			The ratio of the molecular weight of the constituent gas to that assumed
528			for air. Defaults to the ratio for water vapor to dry air
529			(:math:`\epsilon\approx0.622`).
530
531			Returns
532			-------
533			`pint.Quantity`
534			The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)
535
536			Notes
537			-----
538			This function is a straightforward implementation of the equation given in many places,
539			such as [5]_:
540
541			.. math:: r = \epsilon \frac{e}{p - e}
542
543			References
544			----------
545			.. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
546			Survey. 73.
547
548			See Also
549			--------
550			saturation_mixing_ratio, vapor_pressure
551			"""
552			return (molecular_weight_ratio * part_press / (tot_press - part_press)).to("dimensionless")
553
554
555			@exporter.export
556			def saturation_mixing_ratio(tot_press, temperature):
557			r"""Calculate the saturation mixing ratio of water vapor.
558
559			This calculation is given total pressure and the temperature. The implementation
560			uses the formula outlined in [6]_.
561
562			Parameters
563			----------
564			tot_press: `pint.Quantity`
565			Total atmospheric pressure
566			temperature: `pint.Quantity`
567			The temperature
568
569			Returns
570			-------
571			`pint.Quantity`
572			The saturation mixing ratio, dimensionless
573
574			References
575			----------
576			.. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
577			Survey. 73.
578			"""
579			return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)
580
581
582			@exporter.export
583			def equivalent_potential_temperature(pressure, temperature):
584			r"""Calculate equivalent potential temperature.
585
586			This calculation must be given an air parcel's pressure and temperature.
587			The implementation uses the formula outlined in [7]_.
588
589			Parameters
590			----------
591			pressure: `pint.Quantity`
592			Total atmospheric pressure
593			temperature: `pint.Quantity`
594			The temperature
595
596			Returns
597			-------
598			`pint.Quantity`
599			The corresponding equivalent potential temperature of the parcel
600
601			Notes
602			-----
603			.. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}
604
605			References
606			----------
607			.. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
608			Survey. 78-79.
609			"""
610			pottemp = potential_temperature(pressure, temperature)
611			smixr = saturation_mixing_ratio(pressure, temperature)
612			return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))
613

Unidata / MetPy

Pull Request — master (#322)

el() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like