"""Functions to find melting layer from data."""
import numpy as np
import numpy.typing as npt
from atmoslib.constants import T0
from numpy import ma
from scipy.ndimage import gaussian_filter
from cloudnetpy import utils
from cloudnetpy.categorize import droplet
from cloudnetpy.categorize.containers import ClassData
MIN_LAPSE_RATE = 0.004 # K m-1
MIN_LDR_COVERAGE = 0.5
T0_TOLERANCE = 1 # K
[docs]
def find_melting_layer(obs: ClassData, *, smooth: bool = True) -> npt.NDArray:
"""Finds melting layer from model temperature, ldr, and velocity.
Melting layer is detected using linear depolarization ratio, *ldr*,
Doppler velocity, *v*, and wet-bulb temperature, *Tw*.
The algorithm is based on *ldr* having a clear Gaussian peak around
the melting layer. This signature is caused by the growth of ice
crystals into snowflakes that are much larger. In addition, when snow and
ice melt, emerging heavy water droplets start to drop rapidly towards
ground. Thus, there is also a similar positive peak in the
first difference of *v*.
The peak in *ldr* is the primary parameter we analyze. If
*ldr* has a proper peak, and *v* < -1 m/s in the base, melting layer
has been found. If *ldr* is missing we only analyze the behaviour
of *v*, which is always present, to detect the melting layer. If *ldr*
covers the profile but shows no peak, there is no melting layer.
Model temperature is used to limit the melting layer search to a certain
temperature range around 0 C. For ECMWF the range is -4..+3, and for
the rest -8..+6. Model temperature must also reach -1 C in the profile,
and the search is limited in altitude above the highest 0 C level,
because in temperature inversions the temperature range alone may
span several kilometers.
Notes:
This melting layer detection method is novel and needs to be validated.
Also note that there might be some detection problems with strong
updrafts of air. In these cases the absolute values for speed do not
make sense (rain drops can even move upwards instead of down).
Args:
obs: The :class:`ClassData` instance.
smooth: If True, apply a small Gaussian smoother to the
melting layer. Default is True.
Returns:
2-D boolean array denoting the melting layer.
"""
melting_layer = np.zeros(obs.tw.shape, dtype=bool)
ldr_prof: npt.NDArray | None = None
ldr_dprof: npt.NDArray | None = None
ldr_diff: npt.NDArray | None = None
width_prof = None
if hasattr(obs, "ldr"):
# Required for peak detection
diffu = ma.array(np.diff(obs.ldr, axis=1))
ldr_diff = diffu.filled(0)
t_range = _find_model_temperature_range(obs.model_type)
for ind, t_prof in enumerate(obs.tw):
temp_indices = _get_temp_indices(t_prof, t_range, obs.height)
if len(temp_indices) <= 1:
continue
z_prof = obs.z[ind, temp_indices]
v_prof = obs.v[ind, temp_indices]
if ldr_diff is not None:
if not hasattr(obs, "ldr"):
msg = "ldr_diff is not None but obs.ldr does not exist"
raise RuntimeError(msg)
ldr_prof = obs.ldr[ind, temp_indices]
ldr_dprof = ldr_diff[ind, temp_indices]
if (ldr_prof is not None and ma.count(ldr_prof) > 3) or (
v_prof is not None and ma.count(v_prof) > 3
):
try:
if ldr_prof is None or ldr_dprof is None:
msg = "ldr_prof or ldr_dprof is None"
raise AssertionError(msg) # noqa: TRY301
indices = _find_melting_layer_from_ldr(
ldr_prof,
ldr_dprof,
v_prof,
z_prof,
)
except (ValueError, IndexError, AssertionError):
if _has_ldr_coverage(ldr_prof, v_prof):
continue
height = obs.height[temp_indices]
if hasattr(obs, "width"):
width_prof = obs.width[ind, temp_indices]
indices = _find_melting_layer_from_v(v_prof, width_prof, height)
if indices is not None:
melting_layer[ind, temp_indices[indices]] = True
if smooth:
smoothed_layer = gaussian_filter(np.array(melting_layer, dtype=float), (2, 0.1))
melting_layer = (smoothed_layer > 0.2).astype(bool)
return melting_layer
def _has_ldr_coverage(ldr_prof: npt.NDArray | None, v_prof: npt.NDArray) -> bool:
"""Checks if ldr covers enough of the profile to rule out melting layer."""
if ldr_prof is None:
return False
is_v = ~ma.getmaskarray(v_prof)
is_ldr = ~ma.getmaskarray(ldr_prof)
n_ldr = np.count_nonzero(is_ldr & is_v)
return bool(n_ldr >= MIN_LDR_COVERAGE * np.count_nonzero(is_v))
def _find_melting_layer_from_ldr(
ldr_prof: npt.NDArray,
ldr_dprof: npt.NDArray,
v_prof: npt.NDArray,
z_prof: npt.NDArray,
) -> npt.NDArray | None:
peak = int(np.argmax(ldr_prof))
base, top = _basetop(ldr_dprof, peak)
conditions = (
ldr_prof[peak] - ldr_prof[base] > 4,
ldr_prof[peak] > -30,
z_prof[base] > -25,
v_prof[base] < -1.5,
)
if all(conditions):
base = int(np.floor(base + (peak - base) / 2))
return np.arange(base, top)
return None
def _find_melting_layer_from_v(
v_prof: npt.NDArray,
width_prof: npt.NDArray | None,
height: npt.NDArray,
) -> npt.NDArray | None:
v = np.copy(v_prof[:-1])
v_diff = np.diff(v_prof)
v[v_diff < 0] = 0
v[v_diff > 0] = 1
n_increasing = utils.cumsumr(v)
try:
top = int(np.argmax(n_increasing))
base = np.where(n_increasing[:top] == 0)[0][-1]
except IndexError:
return None
if width_prof is not None:
conditions = [
width_prof[base] - width_prof[top] > 0.2,
v_prof[top] - v_prof[base] > 0.5,
50 < (height[top] - height[base]) < 1000,
v_prof[base] < -2,
]
else:
conditions = [
v_prof[top] - v_prof[base] > 2,
50 < (height[top] - height[base]) < 1000,
v_prof[base] < -2,
]
if all(conditions):
base = round(top - (top - base) / 2)
return np.arange(base, top)
return None
def _basetop(dprof: npt.NDArray, pind: int) -> tuple[int, int]:
"""Finds the base and top of ldr peak."""
top = droplet.ind_top(dprof, pind, len(dprof), 10, 2)
base = droplet.ind_base(dprof, pind, 10, 2)
return base, top
def _get_temp_indices(
t_prof: npt.NDArray, t_range: tuple, height: npt.NDArray | None = None
) -> npt.NDArray:
"""Finds indices of temperature profile covering the given range."""
in_range = (t_prof > min(t_range) + T0) & (t_prof < max(t_range) + T0)
if height is not None:
in_range &= height <= _find_max_height(t_prof, t_range, height)
ind = np.where(in_range)[0]
return np.array([]) if len(ind) == 0 else np.arange(np.min(ind), np.max(ind) + 1)
def _find_max_height(t_prof: npt.NDArray, t_range: tuple, height: npt.NDArray) -> float:
"""Finds maximum melting layer height assuming realistic lapse rate."""
warm = np.where(t_prof >= T0)[0]
if len(warm) == 0:
# Model may be slightly too cold when melting layer is near the ground
warm = np.where(t_prof >= T0 - T0_TOLERANCE)[0]
if len(warm) == 0:
return -np.inf
return height[warm[-1]] - min(t_range) / MIN_LAPSE_RATE
def _find_model_temperature_range(model_type: str) -> tuple[float, float]:
"""Returns temperature range around 0C for given model type."""
if model_type == "gdas1":
return -8, 6
return -4, 3