pycontrails.datalib.ecmwf.model_levels

Utilities for working with ECMWF model-level data.

Functions

ml_to_pl(ds, target_pl, *[, lnsp, sp])	Interpolate L137 model-level meteorology data to pressure levels.
model_level_pressure(sp, model_levels)	Return the pressure levels at each model level given the surface pressure.
model_level_reference_pressure([alt_ft_min, ...])	Return the pressure levels at each model level assuming a constant surface pressure.
searchsorted2d(a, v)	Return the indices where elements in v would be inserted in a along its second axis.

pycontrails.datalib.ecmwf.model_levels.ml_to_pl(ds, target_pl, *, lnsp=None, sp=None)

Interpolate L137 model-level meteorology data to pressure levels.

The implementation is here is consistent with ECMWF’s suggested implementation.

Parameters:

• ds (xarray.Dataset) – Dataset with model-level meteorology data. Must include a “model_level” dimension which is not split across chunks. The non-“model_level” dimensions must be aligned with the “lnsp” parameter. Can include any number of variables. Any non-dimension coordinates will be dropped.

• target_pl (npt.ArrayLike) – Target pressure levels, [\(hPa\)].

• lnsp (xarray.DataArray) – Natural logarithm of surface pressure, [\(\ln(\text{Pa})\)]. If provided, sp is ignored. At least one of lnsp or sp must be provided. The chunking over dimensions in common with ds must be the same as ds.

• sp (xarray.DataArray) – Surface pressure, [\(\text{Pa}\)]. At least one of lnsp or sp must be provided. The chunking over dimensions in common with ds must be the same as ds.

Returns:

xarray.Dataset – Interpolated data on the target pressure levels. This has the same dimensions as ds except that the “model_level” dimension is replaced with “level”. The shape of the “level” dimension is the length of target_pl. If ds is dask-backed, the output will be as well. Call .compute() to compute the result eagerly.

pycontrails.datalib.ecmwf.model_levels.model_level_pressure(sp, model_levels)

Return the pressure levels at each model level given the surface pressure.

This function assumes 137 model levels. Unlike model_level_reference_pressure(), this function does not assume constant pressure. Instead, it uses the half-level pressure formula \(p = a + b \cdot \text{sp}\) where \(a\) and \(b\) are constants for each model level.

Parameters:

• sp (xarray.DataArray) – Surface pressure, [\(\text{Pa}\)]. A warning is issued if the minimum value of sp is less than 30320.0 Pa. Such low values are unrealistic.

• model_levels (npt.ArrayLike) – Target model levels. Expected to be a one-dimensional array of integers between 1 and 137.

Returns:

xarray.DataArray – Pressure levels at each model level, [\(hPa\)]. The shape of the output is the product of the shape of the input and the length of model_levels. In other words, the output will have dimensions of the input plus a new dimension for model_levels.

If sp is not dask-backed, the output will be computed eagerly. In particular, if sp has a large size and model_levels is a large range, this function may consume a large amount of memory.

The dtype of the output is the same as the dtype of the sp parameter.

Examples

>>> import numpy as np
>>> import xarray as xr

>>> sp_arr = np.linspace(101325.0, 90000.0, 16).reshape(4, 4)
>>> longitude = np.linspace(-180, 180, 4)
>>> latitude = np.linspace(-90, 90, 4)
>>> sp = xr.DataArray(sp_arr, coords={"longitude": longitude, "latitude": latitude})

>>> model_levels = [80, 100]
>>> model_level_pressure(sp, model_levels)
<xarray.DataArray (model_level: 2, longitude: 4, latitude: 4)> Size: 256B
array([[[259.75493944, 259.27107504, 258.78721064, 258.30334624],
        [257.81948184, 257.33561744, 256.85175304, 256.36788864],
        [255.88402424, 255.40015984, 254.91629544, 254.43243104],
        [253.94856664, 253.46470224, 252.98083784, 252.49697344]],
       [[589.67975444, 586.47283154, 583.26590864, 580.05898574],
        [576.85206284, 573.64513994, 570.43821704, 567.23129414],
        [564.02437124, 560.81744834, 557.61052544, 554.40360254],
        [551.19667964, 547.98975674, 544.78283384, 541.57591094]]])
Coordinates:
  * longitude    (longitude) float64 32B -180.0 -60.0 60.0 180.0
  * latitude     (latitude) float64 32B -90.0 -30.0 30.0 90.0
  * model_level  (model_level) int64 16B 80 100

See also

model_level_reference_pressure

pycontrails.datalib.ecmwf.model_levels.model_level_reference_pressure(alt_ft_min=None, alt_ft_max=None)

Return the pressure levels at each model level assuming a constant surface pressure.

This function assumes 137 model levels and the constant ICAO ISA surface pressure of 1013.25 hPa.

The returned pressure levels are rounded to the nearest hPa.

Parameters:

• alt_ft_min (float | None) – Minimum altitude, [\(ft\)]. If None, there is no minimum altitude used in filtering the MODEL_LEVELS_PATH table.

• alt_ft_max (float | None) – Maximum altitude, [\(ft\)]. If None, there is no maximum altitude used in filtering the MODEL_LEVELS_PATH table.

Returns:

list[int] – List of pressure levels, [\(hPa\)] between the minimum and maximum altitudes.

See also

model_level_pressure

pycontrails.datalib.ecmwf.model_levels.searchsorted2d(a, v)

Return the indices where elements in v would be inserted in a along its second axis.

Implementation based on a StackOverflow answer.

Parameters:

• a (npt.NDArray[np.floating]) – 2D array of shape (m, n) that is sorted along its second axis. This is not checked.

• v (npt.NDArray[np.floating]) – 1D array of values of shape (k,) to insert into the second axis of a. The current implementation could be extended to handle 2D arrays as well.

Returns:

npt.NDArray[np.int64] – 2D array of indices where elements in v would be inserted in a along its second axis to keep the second axis of a sorted. The shape of the output is (m, k).

Examples

>>> a = np.array([
...  [ 1.,  8., 11., 12.],
...  [ 5.,  8.,  9., 14.],
...  [ 4.,  5.,  6., 17.],
...  ])
>>> v = np.array([3., 7., 10., 13., 15.])
>>> searchsorted2d(a, v)
array([[1, 1, 2, 4, 4],
       [0, 1, 3, 3, 4],
       [0, 3, 3, 3, 3]])