cape() - 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 23:39 UTC

cape() A

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	2
dl	0
loc	60
rs	9.5555
c	1
b	0
f	0

How to fix Long Method

# 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):
    r"""Calculate convective available potential energy (CAPE).

        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
        temperature_parcel : `pint.Quantity`
            The temperature of the parcel

        Returns
        -------
        `pint.Quantity`
            Convective available potential energy (CAPE).

        See Also
        --------
        virtual_temperature, find_intersections, lfc, el
        """
    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			r"""Calculate convective available potential energy (CAPE).
342
343			Description
344
345			Parameters
346			----------
347			pressure : `pint.Quantity`
348			The atmospheric pressure level(s) of interest. The first entry should be the starting
349			point pressure.
350			temperature : `pint.Quantity`
351			The starting temperature
352			dewpt : `pint.Quantity`
353			The starting dew point
354			temperature_parcel : `pint.Quantity`
355			The temperature of the parcel
356
357			Returns
358			-------
359			`pint.Quantity`
360			Convective available potential energy (CAPE).
361
362			See Also
363			--------
364			virtual_temperature, find_intersections, lfc, el
365			"""
366			if virtual_temp:
367			temperature = virtual_temperature(pressure, temperature, dewpt)
368			temperature_parcel = virtual_temperature(pressure, temperature_parcel, dewpt)
369
370			x = np.copy(pressure)
371			y = temperature_parcel - temperature
372
373			# Find and append crossings to the data
374			crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
375			x = np.concatenate((x, crossings[0]))
376			y = np.concatenate((y, crossings[1]))
377
378			# Resort so that data is in pressure order
379			sort_idx = np.argsort(x)
380			x = x[sort_idx]
381			y = y[sort_idx]
382
383			# Remove duplicate data points if there are any
384			keep_idx = np.ediff1d(x, to_end=[1]) >= 1
385			x = x[keep_idx]
386			y = y[keep_idx]
387
388			# Clip out negative area (temperature parcel < temperature environment)
389			y[y <= 0] = 0
390
391			# Only use data between the LFC and EL for calculation
392			lfc_pressure = lfc(pressure, temperature, dewpt)[0].magnitude
393			el_pressure = el(pressure, temperature, dewpt)[0].magnitude
394			p_mask = (x <= lfc_pressure) & (x >= el_pressure)
395			x = x[p_mask]
396			y = y[p_mask]
397
398			return Rd * np.trapz(y, np.log(x))
399
400
401			@exporter.export
402			def vapor_pressure(pressure, mixing):
403			r"""Calculate water vapor (partial) pressure.
404
405			Given total `pressure` and water vapor `mixing` ratio, calculates the
406			partial pressure of water vapor.
407
408			Parameters
409			----------
410			pressure : `pint.Quantity`
411			total atmospheric pressure
412			mixing : `pint.Quantity`
413			dimensionless mass mixing ratio
414
415			Returns
416			-------
417			`pint.Quantity`
418			The ambient water vapor (partial) pressure in the same units as
419			`pressure`.
420
421			Notes
422			-----
423			This function is a straightforward implementation of the equation given in many places,
424			such as [2]_:
425
426			.. math:: e = p \frac{r}{r + \epsilon}
427
428			References
429			----------
430			.. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
431			Survey. 71.
432
433			See Also
434			--------
435			saturation_vapor_pressure, dewpoint
436			"""
437			return pressure * mixing / (epsilon + mixing)
438
439
440			@exporter.export
441			def saturation_vapor_pressure(temperature):
442			r"""Calculate the saturation water vapor (partial) pressure.
443
444			Parameters
445			----------
446			temperature : `pint.Quantity`
447			The temperature
448
449			Returns
450			-------
451			`pint.Quantity`
452			The saturation water vapor (partial) pressure
453
454			See Also
455			--------
456			vapor_pressure, dewpoint
457
458			Notes
459			-----
460			Instead of temperature, dewpoint may be used in order to calculate
461			the actual (ambient) water vapor (partial) pressure.
462
463			The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:
464
465			.. math:: 6.112 e^\frac{17.67T}{T + 243.5}
466
467			References
468			----------
469			.. [3] Bolton, D., 1980: The Computation of Equivalent Potential
470			Temperature. Mon. Wea. Rev., 108, 1046-1053.
471			"""
472			# Converted from original in terms of C to use kelvin. Using raw absolute values of C in
473			# a formula plays havoc with units support.
474			return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
475			(temperature - 29.65 * units.kelvin))
476
477
478			@exporter.export
479			def dewpoint_rh(temperature, rh):
480			r"""Calculate the ambient dewpoint given air temperature and relative humidity.
481
482			Parameters
483			----------
484			temperature : `pint.Quantity`
485			Air temperature
486			rh : `pint.Quantity`
487			Relative humidity expressed as a ratio in the range [0, 1]
488
489			Returns
490			-------
491			`pint.Quantity`
492			The dew point temperature
493
494			See Also
495			--------
496			dewpoint, saturation_vapor_pressure
497			"""
498			return dewpoint(rh * saturation_vapor_pressure(temperature))
499
500
501			@exporter.export
502			def dewpoint(e):
503			r"""Calculate the ambient dewpoint given the vapor pressure.
504
505			Parameters
506			----------
507			e : `pint.Quantity`
508			Water vapor partial pressure
509
510			Returns
511			-------
512			`pint.Quantity`
513			Dew point temperature
514
515			See Also
516			--------
517			dewpoint_rh, saturation_vapor_pressure, vapor_pressure
518
519			Notes
520			-----
521			This function inverts the Bolton 1980 [4]_ formula for saturation vapor
522			pressure to instead calculate the temperature. This yield the following
523			formula for dewpoint in degrees Celsius:
524
525			.. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}
526
527			References
528			----------
529			.. [4] Bolton, D., 1980: The Computation of Equivalent Potential
530			Temperature. Mon. Wea. Rev., 108, 1046-1053.
531			"""
532			val = np.log(e / sat_pressure_0c)
533			return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
534
535
536			@exporter.export
537			def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
538			r"""Calculate the mixing ratio of a gas.
539
540			This calculates mixing ratio given its partial pressure and the total pressure of
541			the air. There are no required units for the input arrays, other than that
542			they have the same units.
543
544			Parameters
545			----------
546			part_press : `pint.Quantity`
547			Partial pressure of the constituent gas
548			tot_press : `pint.Quantity`
549			Total air pressure
550			molecular_weight_ratio : `pint.Quantity` or float, optional
551			The ratio of the molecular weight of the constituent gas to that assumed
552			for air. Defaults to the ratio for water vapor to dry air
553			(:math:`\epsilon\approx0.622`).
554
555			Returns
556			-------
557			`pint.Quantity`
558			The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)
559
560			Notes
561			-----
562			This function is a straightforward implementation of the equation given in many places,
563			such as [5]_:
564
565			.. math:: r = \epsilon \frac{e}{p - e}
566
567			References
568			----------
569			.. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
570			Survey. 73.
571
572			See Also
573			--------
574			saturation_mixing_ratio, vapor_pressure
575			"""
576			return (molecular_weight_ratio * part_press / (tot_press - part_press)).to("dimensionless")
577
578
579			@exporter.export
580			def saturation_mixing_ratio(tot_press, temperature):
581			r"""Calculate the saturation mixing ratio of water vapor.
582
583			This calculation is given total pressure and the temperature. The implementation
584			uses the formula outlined in [6]_.
585
586			Parameters
587			----------
588			tot_press: `pint.Quantity`
589			Total atmospheric pressure
590			temperature: `pint.Quantity`
591			The temperature
592
593			Returns
594			-------
595			`pint.Quantity`
596			The saturation mixing ratio, dimensionless
597
598			References
599			----------
600			.. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
601			Survey. 73.
602			"""
603			return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)
604
605
606			@exporter.export
607			def equivalent_potential_temperature(pressure, temperature):
608			r"""Calculate equivalent potential temperature.
609
610			This calculation must be given an air parcel's pressure and temperature.
611			The implementation uses the formula outlined in [7]_.
612
613			Parameters
614			----------
615			pressure: `pint.Quantity`
616			Total atmospheric pressure
617			temperature: `pint.Quantity`
618			The temperature
619
620			Returns
621			-------
622			`pint.Quantity`
623			The corresponding equivalent potential temperature of the parcel
624
625			Notes
626			-----
627			.. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}
628
629			References
630			----------
631			.. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
632			Survey. 78-79.
633			"""
634			pottemp = potential_temperature(pressure, temperature)
635			smixr = saturation_mixing_ratio(pressure, temperature)
636			return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))
637

Unidata / MetPy

Pull Request — master (#322)

cape() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like