add lifted_index and precipitable_water functions in meteolib package

meteoinfo-lab/milconfig.xml

@@ -13,18 +13,18 @@
     <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
     <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\interpolate"/>
     <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array\slice"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
     <RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
     <RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\calc"/>
   </Path>
   <File>
     <OpenedFiles>
-      <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\mixed_layer_cape_cin.py"/>
-      <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\isentropic_analysis.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\lifted_index.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\precipitable_water.py"/>
     </OpenedFiles>
     <RecentFiles>
-      <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\mixed_layer_cape_cin.py"/>
-      <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\isentropic_analysis.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\lifted_index.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\precipitable_water.py"/>
     </RecentFiles>
   </File>
   <Font>

meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__.py

@@ -2,9 +2,11 @@ from .kinematics import *
 from .thermo import *
 from .tools import *
 from .basic import *
+from .indices import *
 
 __all__ = []
 __all__ += kinematics.__all__
 __all__ += thermo.__all__
 __all__ += tools.__all__
-__all__ += basic.__all__
+__all__ += basic.__all__
+__all__ += indices.__all__

meteoinfo-lab/pylib/mipylib/meteolib/calc/indices.py

@@ -0,0 +1,89 @@
+"""
+Contains calculation of various derived indices.
+
+Ported from MetPy.
+"""
+
+import mipylib.numeric as np
+from .tools import _remove_nans, get_layer
+from .thermo import mixing_ratio, saturation_vapor_pressure
+from .. import constants
+
+
+__all__ = ['precipitable_water']
+
+
+def precipitable_water(pressure, dewpoint, bottom=None, top=None):
+    r"""Calculate precipitable water through the depth of a sounding.
+
+    Formula used is:
+
+    .. math::  -\frac{1}{\rho_l g} \int\limits_{p_\text{bottom}}^{p_\text{top}} r dp
+
+    from [Salby1996]_, p. 28
+
+    Parameters
+    ----------
+    pressure : `hPa`
+        Atmospheric pressure profile
+
+    dewpoint : `kelvin`
+        Atmospheric dewpoint profile
+
+    bottom: `hPa`, optional
+        Bottom of the layer, specified in pressure. Defaults to None (highest pressure).
+
+    top: `hPa`, optional
+        Top of the layer, specified in pressure. Defaults to None (lowest pressure).
+
+    Returns
+    -------
+    `millimeters`
+        Precipitable water in the layer
+
+    Examples
+    --------
+    >>> from mipylib.meteolib.calc import precipitable_water
+    >>> from mipylib.meteolib import constants as cons
+    >>> pressure = np.array([1000, 950, 900])
+    >>> dewpoint = np.array([20, 15, 10]) + cons.degCtoK
+    >>> pw = precipitable_water(pressure, dewpoint)
+    (11.7702606153 millimeter)
+
+    Notes
+    -----
+    Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
+
+    """
+    # Sort pressure and dewpoint to be in decreasing pressure order (increasing height)
+    sort_inds = np.argsort(pressure)[::-1]
+    pressure = pressure[sort_inds]
+    dewpoint = dewpoint[sort_inds]
+
+    pressure, dewpoint = _remove_nans(pressure, dewpoint)
+
+    min_pressure = np.min(pressure)
+    max_pressure = np.max(pressure)
+
+    if top is None:
+        top = min_pressure
+    elif not min_pressure <= top <= max_pressure:
+        raise ValueError('The pressure and dewpoint profile ranges from {} to '
+                         '{}, after removing missing values. {} is outside '
+                         'this range.'.format(max_pressure, min_pressure, top))
+
+    if bottom is None:
+        bottom = max_pressure
+    elif not min_pressure <= bottom <= max_pressure:
+        raise ValueError('The pressure and dewpoint profile ranges from {} to '
+                         '{}, after removing missing values. {} is outside '
+                         'this range.'.format(max_pressure, min_pressure, bottom))
+
+    pres_layer, dewpoint_layer = get_layer(pressure, dewpoint, bottom=bottom,
+                                           depth=bottom - top)
+
+    w = mixing_ratio(saturation_vapor_pressure(dewpoint_layer), pres_layer)
+
+    # Since pressure is in decreasing order, pw will be the opposite sign of that expected.
+    pw = -np.trapezoid(w, pres_layer) / (constants.g * constants.rho_l)
+    return pw * 1e5

meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo.py

@@ -20,8 +20,8 @@ from mipylib.numeric.special import lambertw
 __all__ = [
     'cape_cin','ccl','density','dewpoint','dewpoint_from_relative_humidity','dry_lapse','dry_static_energy',
     'el','equivalent_potential_temperature','exner_function',
-    'isentropic_interpolation','lcl','lfc','mixed_layer','mixed_layer_cape_cin','mixed_parcel','mixing_ratio',
-    'mixing_ratio_from_relative_humidity',
+    'isentropic_interpolation','lcl','lfc','lifted_index','mixed_layer','mixed_layer_cape_cin',
+    'mixed_parcel','mixing_ratio','mixing_ratio_from_relative_humidity',
     'mixing_ratio_from_specific_humidity',
     'moist_air_gas_constant','moist_air_specific_heat_pressure',
     'moist_air_poisson_exponent','moist_lapse','most_unstable_parcel','most_unstable_cape_cin',
@@ -838,6 +838,8 @@ def exner_function(pressure, reference_pressure=constants.P0):
     potential_temperature
     temperature_from_potential_temperature
     """
+    if isinstance(pressure, (list, tuple)):
+        pressure = np.asarray(pressure)
     return (pressure / reference_pressure)**constants.kappa
 
 
@@ -3070,6 +3072,87 @@ def specific_humidity_from_relative_humidity(pressure, temperature, rh):
     return specific_humidity_from_dewpoint(pressure, dp)
 
 
+def lifted_index(pressure, temperature, parcel_profile, vertical_dim=0):
+    """Calculate Lifted Index from the pressure temperature and parcel profile.
+
+    Lifted index formula derived from [Galway1956]_ and referenced by [DoswellSchultz2006]_:
+
+    .. math:: LI = T500 - Tp500
+
+    where:
+
+    * :math:`T500` is the measured temperature at 500 hPa
+    * :math:`Tp500` is the temperature of the lifted parcel at 500 hPa
+
+    Calculation of the lifted index is defined as the temperature difference between the
+    observed 500 hPa temperature and the temperature of a parcel lifted from the
+    surface to 500 hPa. As noted in [Galway1956]_, a low-level mixed parcel is often used
+    as the surface parcel.
+
+    Parameters
+    ----------
+    pressure : `hPa`
+        Atmospheric pressure level(s) of interest, in order from highest to
+        lowest pressure
+
+    temperature : `kelvin`
+        Atmospheric temperature corresponding to pressure
+
+    parcel_profile : `kelvin`
+        Temperature profile of the parcel
+
+    vertical_dim : int, optional
+        The axis corresponding to vertical, defaults to 0. Automatically determined from
+        xarray DataArray arguments.
+
+    Returns
+    -------
+    `delta temperature`
+        Lifted Index
+
+    Examples
+    --------
+    >>> from mipylib.meteolib.calc import lifted_index, mixed_parcel, parcel_profile
+    >>> from mipylib.meteolib import constants as cons
+    >>>
+    >>> # Define pressure, temperature, dewpoint temperature, and height
+    >>> p = [1002., 1000., 993., 925., 850., 846., 723., 632., 479., 284.,
+    ...      239., 200., 131., 91., 72.7, 54.6, 41., 30., 22.8]
+    >>> T = array([28.2, 27., 25.4, 19.4, 12.8, 12.3, 4.2, 0.8, -12.7, -41.7, -52.3,
+    ...      -57.5, -54.9, -57.8, -58.5, -52.3, -53.4, -50.3, -49.9]) + cons.degCtoK
+    >>> Td = array([14.2, 14., 12.4, 11.4, 10.2, 10.1, -7.8, -16.2, -37.7, -55.7,
+    ...       -58.3, -69.5, -85.5, -88., -89.5, -88.3, -88.4, -87.3, -87.9]) + cons.degCtoK
+    >>> h = [139, 159, 221, 839, 1559, 1599, 2895, 3982, 6150, 9933, 11072,
+    ...      12200, 14906, 17231, 18650, 20474, 22323, 24350, 26149]
+    >>>
+    >>> # Calculate 500m mixed parcel
+    >>> parcel_p, parcel_t, parcel_td = mixed_parcel(p, T, Td, depth=500, height=h)
+    >>>
+    >>> # Replace sounding temp/pressure in lowest 500m with mixed values
+    >>> above = h > 500
+    >>> press = np.concatenate([[parcel_p], p[above]])
+    >>> temp = np.concatenate([[parcel_t], T[above]])
+    >>>
+    >>> # Calculate parcel profile from our new mixed parcel
+    >>> mixed_prof = parcel_profile(press, parcel_t, parcel_td)
+    >>>
+    >>> # Calculate lifted index using our mixed profile
+    >>> lifted_index(press, temp, mixed_prof)
+    <([2.54198585], 'delta_degree_Celsius')>
+
+    See Also
+    --------
+    mixed_parcel, parcel_profile
+
+    """
+    # find the measured temperature and parcel profile temperature at 500 hPa.
+    t500, tp500 = interpolate_1d(500,
+                                 pressure, temperature, parcel_profile, axis=vertical_dim)
+
+    # calculate the lifted index.
+    return t500 - tp500
+
+
 def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
                  height=None, bottom=None, depth=None, interpolate=True):
     r"""Calculate the properties of a parcel mixed from a layer.
@@ -3146,6 +3229,10 @@ def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
     Quantities even when given xarray DataArray profiles.
 
     """
+    pressure = np.asarray(pressure)
+    temperature = np.asarray(temperature)
+    dewpoint = np.asarray(dewpoint)
+
     # If a parcel starting pressure is not provided, use the surface
     if not parcel_start_pressure:
         parcel_start_pressure = pressure[0]
@@ -3230,7 +3317,7 @@ def mixed_layer(pressure, *args, **kwargs):
     >>> # calculate dewpoint
     >>> Td = dewpoint_from_relative_humidity(T, rh)
     >>> # find mixed layer T and Td of depth 50 hPa
-    >>> mixed_layer(p, T, Td, depth=50 * units.hPa)
+    >>> mixed_layer(p, T, Td, depth=50)
     [<(26.5798571, 'degree_Celsius')>, <(16.6455209, 'degree_Celsius')>]
 
     Notes

meteoinfo-lab/pylib/mipylib/meteolib/calc/tools.py

@@ -10,8 +10,9 @@ from mipylib.geolib import Geod
 from ..interpolate import interpolate_1d, log_interpolate_1d
 from ..cbook import broadcast_indices
 
-__all__ = ['resample_nn_1d', 'nearest_intersection_idx', 'first_derivative', 'find_intersections', 'gradient',
-           'lat_lon_grid_deltas', 'get_layer_heights', 'find_bounding_indices', 'geospatial_gradient']
+__all__ = ['resample_nn_1d', 'nearest_intersection_idx', 'first_derivative', 'find_intersections', 'get_layer',
+           'gradient', 'lat_lon_grid_deltas', 'get_layer_heights', 'find_bounding_indices',
+           'geospatial_gradient']
 
 
 def resample_nn_1d(a, centers):
@@ -276,8 +277,8 @@ def _get_bound_pressure_height(pressure, bound, height=None, interpolate=True, i
                     # Need to cast back to the input type since interp (up to at least numpy
                     # 1.13 always returns float64. This can cause upstream users problems,
                     # resulting in something like np.append() to upcast.
-                    bound_pressure = np.interp(np.atleast_1d(bound),
-                                               height, pressure).astype(np.result_type(bound))
+                    bound_pressure = interpolate_1d(np.atleast_1d(bound),
+                                               height, pressure).astype(bound.dtype)
                 else:
                     idx = (np.abs(height - bound)).argmin()
                     bound_pressure = pressure[idx]
@@ -435,7 +436,7 @@ def get_layer(pressure, *args, **kwargs):
     bottom = kwargs.pop('bottom', None)
     depth = kwargs.pop('depth', 100)    #'hPa'
     interpolate = kwargs.pop('interpolate', True)
-    is_pressure = kwargs.pop('is_pressure', True)
+    is_pressure = height is None
 
     # Make sure pressure and datavars are the same length
     for datavar in args:
@@ -457,7 +458,7 @@ def get_layer(pressure, *args, **kwargs):
         top = bottom_height + depth
 
     top_pressure, _ = _get_bound_pressure_height(pressure, top, height=height,
-                                                 interpolate=interpolate)
+                                                 interpolate=interpolate, is_pressure=is_pressure)
 
     ret = []  # returned data variables in layer