diff --git a/meteoinfo-lab/milconfig.xml b/meteoinfo-lab/milconfig.xml
index 99beaeab..fcd9b3e3 100644
--- a/meteoinfo-lab/milconfig.xml
+++ b/meteoinfo-lab/milconfig.xml
@@ -13,18 +13,18 @@
-
+
-
-
+
+
-
-
+
+
diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__$py.class b/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__$py.class
index 9a6484c8..894615d5 100644
Binary files a/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__$py.class and b/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__$py.class differ
diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__.py b/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__.py
index 5b8d8425..db067643 100644
--- a/meteoinfo-lab/pylib/mipylib/meteolib/calc/__init__.py
+++ b/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__
\ No newline at end of file
+__all__ += basic.__all__
+__all__ += indices.__all__
diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/calc/indices.py b/meteoinfo-lab/pylib/mipylib/meteolib/calc/indices.py
new file mode 100644
index 00000000..3c0eac4d
--- /dev/null
+++ b/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
diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo$py.class b/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo$py.class
index e07d6456..cf0c39a1 100644
Binary files a/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo$py.class and b/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo$py.class differ
diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo.py b/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo.py
index 7e9235d3..2ca9cc6b 100644
--- a/meteoinfo-lab/pylib/mipylib/meteolib/calc/thermo.py
+++ b/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
diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/calc/tools.py b/meteoinfo-lab/pylib/mipylib/meteolib/calc/tools.py
index 60462e25..8561bb9e 100644
--- a/meteoinfo-lab/pylib/mipylib/meteolib/calc/tools.py
+++ b/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