virtual_temperature() - 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-16 00:46 UTC

virtual_temperature() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	1
c	1
b	0
f	0
dl	0
loc	27
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):
    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).to(units.degK)

    # Find and append crossings to the data
    crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
    x = concatenate((x, crossings[0]))
    y = 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 * units.degK] = 0 * units.degK

    # 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)) * units.degK)).to(units('J/kg'))


@exporter.export

def cin(pressure, temperature, dewpt, temperature_parcel, virtual_temp=False):
    r"""Calculate convective inhibition (CIN).

        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 inhibition (CIN).

        See Also
        --------
        virtual_temperature, find_intersections, lfc
        """
    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).to(units.degK)

    # Find and append crossings to the data
    crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
    x = concatenate((x, crossings[0]))
    y = 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 positive area (temperature parcel > temperature environment)
    y[y >= 0 * units.degK] = 0 * units.degK

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

    return (Rd * (np.trapz(y, np.log(x)) * units.degK)).to(units('J/kg'))


@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	View Code Duplication	@exporter.export
		0 ignored issues – show Duplication introduced 2017-02-16 00:47 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
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).to(units.degK)
372
373		# Find and append crossings to the data
374		crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
375		x = concatenate((x, crossings[0]))
376		y = 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 * units.degK] = 0 * units.degK
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)) * units.degK)).to(units('J/kg'))
399
400
401	View Code Duplication	@exporter.export
		0 ignored issues – show Duplication introduced 2017-02-16 00:47 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
402		def cin(pressure, temperature, dewpt, temperature_parcel, virtual_temp=False):
403		r"""Calculate convective inhibition (CIN).
404
405		Description
406
407		Parameters
408		----------
409		pressure : `pint.Quantity`
410		The atmospheric pressure level(s) of interest. The first entry should be the starting
411		point pressure.
412		temperature : `pint.Quantity`
413		The starting temperature
414		dewpt : `pint.Quantity`
415		The starting dew point
416		temperature_parcel : `pint.Quantity`
417		The temperature of the parcel
418
419		Returns
420		-------
421		`pint.Quantity`
422		Convective inhibition (CIN).
423
424		See Also
425		--------
426		virtual_temperature, find_intersections, lfc
427		"""
428		if virtual_temp:
429		temperature = virtual_temperature(pressure, temperature, dewpt)
430		temperature_parcel = virtual_temperature(pressure, temperature_parcel, dewpt)
431
432		x = np.copy(pressure)
433		y = (temperature_parcel - temperature).to(units.degK)
434
435		# Find and append crossings to the data
436		crossings = find_intersections(x[1:], y[1:], np.zeros_like(y[1:]))
437		x = concatenate((x, crossings[0]))
438		y = concatenate((y, crossings[1]))
439
440		# Resort so that data is in pressure order
441		sort_idx = np.argsort(x)
442		x = x[sort_idx]
443		y = y[sort_idx]
444
445		# Remove duplicate data points if there are any
446		keep_idx = np.ediff1d(x, to_end=[1]) >= 1
447		x = x[keep_idx]
448		y = y[keep_idx]
449
450		# Clip out positive area (temperature parcel > temperature environment)
451		y[y >= 0 * units.degK] = 0 * units.degK
452
453		# Only use data between the surface and the LFC for calculation
454		lfc_pressure = lfc(pressure, temperature, dewpt)[0].magnitude
455		p_mask = (x >= lfc_pressure)
456		x = x[p_mask]
457		y = y[p_mask]
458
459		return (Rd * (np.trapz(y, np.log(x)) * units.degK)).to(units('J/kg'))
460
461
462		@exporter.export
463		def vapor_pressure(pressure, mixing):
464		r"""Calculate water vapor (partial) pressure.
465
466		Given total `pressure` and water vapor `mixing` ratio, calculates the
467		partial pressure of water vapor.
468
469		Parameters
470		----------
471		pressure : `pint.Quantity`
472		total atmospheric pressure
473		mixing : `pint.Quantity`
474		dimensionless mass mixing ratio
475
476		Returns
477		-------
478		`pint.Quantity`
479		The ambient water vapor (partial) pressure in the same units as
480		`pressure`.
481
482		Notes
483		-----
484		This function is a straightforward implementation of the equation given in many places,
485		such as [2]_:
486
487		.. math:: e = p \frac{r}{r + \epsilon}
488
489		References
490		----------
491		.. [2] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
492		Survey. 71.
493
494		See Also
495		--------
496		saturation_vapor_pressure, dewpoint
497		"""
498		return pressure * mixing / (epsilon + mixing)
499
500
501		@exporter.export
502		def saturation_vapor_pressure(temperature):
503		r"""Calculate the saturation water vapor (partial) pressure.
504
505		Parameters
506		----------
507		temperature : `pint.Quantity`
508		The temperature
509
510		Returns
511		-------
512		`pint.Quantity`
513		The saturation water vapor (partial) pressure
514
515		See Also
516		--------
517		vapor_pressure, dewpoint
518
519		Notes
520		-----
521		Instead of temperature, dewpoint may be used in order to calculate
522		the actual (ambient) water vapor (partial) pressure.
523
524		The formula used is that from Bolton 1980 [3]_ for T in degrees Celsius:
525
526		.. math:: 6.112 e^\frac{17.67T}{T + 243.5}
527
528		References
529		----------
530		.. [3] Bolton, D., 1980: The Computation of Equivalent Potential
531		Temperature. Mon. Wea. Rev., 108, 1046-1053.
532		"""
533		# Converted from original in terms of C to use kelvin. Using raw absolute values of C in
534		# a formula plays havoc with units support.
535		return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) /
536		(temperature - 29.65 * units.kelvin))
537
538
539		@exporter.export
540		def dewpoint_rh(temperature, rh):
541		r"""Calculate the ambient dewpoint given air temperature and relative humidity.
542
543		Parameters
544		----------
545		temperature : `pint.Quantity`
546		Air temperature
547		rh : `pint.Quantity`
548		Relative humidity expressed as a ratio in the range [0, 1]
549
550		Returns
551		-------
552		`pint.Quantity`
553		The dew point temperature
554
555		See Also
556		--------
557		dewpoint, saturation_vapor_pressure
558		"""
559		return dewpoint(rh * saturation_vapor_pressure(temperature))
560
561
562		@exporter.export
563		def dewpoint(e):
564		r"""Calculate the ambient dewpoint given the vapor pressure.
565
566		Parameters
567		----------
568		e : `pint.Quantity`
569		Water vapor partial pressure
570
571		Returns
572		-------
573		`pint.Quantity`
574		Dew point temperature
575
576		See Also
577		--------
578		dewpoint_rh, saturation_vapor_pressure, vapor_pressure
579
580		Notes
581		-----
582		This function inverts the Bolton 1980 [4]_ formula for saturation vapor
583		pressure to instead calculate the temperature. This yield the following
584		formula for dewpoint in degrees Celsius:
585
586		.. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}
587
588		References
589		----------
590		.. [4] Bolton, D., 1980: The Computation of Equivalent Potential
591		Temperature. Mon. Wea. Rev., 108, 1046-1053.
592		"""
593		val = np.log(e / sat_pressure_0c)
594		return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
595
596
597		@exporter.export
598		def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon):
599		r"""Calculate the mixing ratio of a gas.
600
601		This calculates mixing ratio given its partial pressure and the total pressure of
602		the air. There are no required units for the input arrays, other than that
603		they have the same units.
604
605		Parameters
606		----------
607		part_press : `pint.Quantity`
608		Partial pressure of the constituent gas
609		tot_press : `pint.Quantity`
610		Total air pressure
611		molecular_weight_ratio : `pint.Quantity` or float, optional
612		The ratio of the molecular weight of the constituent gas to that assumed
613		for air. Defaults to the ratio for water vapor to dry air
614		(:math:`\epsilon\approx0.622`).
615
616		Returns
617		-------
618		`pint.Quantity`
619		The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g)
620
621		Notes
622		-----
623		This function is a straightforward implementation of the equation given in many places,
624		such as [5]_:
625
626		.. math:: r = \epsilon \frac{e}{p - e}
627
628		References
629		----------
630		.. [5] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
631		Survey. 73.
632
633		See Also
634		--------
635		saturation_mixing_ratio, vapor_pressure
636		"""
637		return (molecular_weight_ratio * part_press / (tot_press - part_press)).to("dimensionless")
638
639
640		@exporter.export
641		def saturation_mixing_ratio(tot_press, temperature):
642		r"""Calculate the saturation mixing ratio of water vapor.
643
644		This calculation is given total pressure and the temperature. The implementation
645		uses the formula outlined in [6]_.
646
647		Parameters
648		----------
649		tot_press: `pint.Quantity`
650		Total atmospheric pressure
651		temperature: `pint.Quantity`
652		The temperature
653
654		Returns
655		-------
656		`pint.Quantity`
657		The saturation mixing ratio, dimensionless
658
659		References
660		----------
661		.. [6] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
662		Survey. 73.
663		"""
664		return mixing_ratio(saturation_vapor_pressure(temperature), tot_press)
665
666
667		@exporter.export
668		def equivalent_potential_temperature(pressure, temperature):
669		r"""Calculate equivalent potential temperature.
670
671		This calculation must be given an air parcel's pressure and temperature.
672		The implementation uses the formula outlined in [7]_.
673
674		Parameters
675		----------
676		pressure: `pint.Quantity`
677		Total atmospheric pressure
678		temperature: `pint.Quantity`
679		The temperature
680
681		Returns
682		-------
683		`pint.Quantity`
684		The corresponding equivalent potential temperature of the parcel
685
686		Notes
687		-----
688		.. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T}
689
690		References
691		----------
692		.. [7] Hobbs, Peter V. and Wallace, John M., 1977: Atmospheric Science, an Introductory
693		Survey. 78-79.
694		"""
695		pottemp = potential_temperature(pressure, temperature)
696		smixr = saturation_mixing_ratio(pressure, temperature)
697		return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))
698

Unidata / MetPy

Pull Request — master (#322)

virtual_temperature() B

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like