"""Module to find insects from data."""
import numpy as np
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
[docs]
def find_insects(
obs: ClassData,
melting_layer: np.ndarray,
liquid_layers: np.ndarray,
prob_lim: float = 0.8,
) -> tuple[np.ndarray, np.ndarray]:
"""Returns insect probability and boolean array of insect presence.
Insects are classified by estimating heuristic probability
of insects from various individual radar parameters and combining
these probabilities. Insects typically yield small echo and spectral width
but high linear depolarization ratio (ldr), and they are present in warm
temperatures.
The combination of echo, ldr and temperature is generally the best proxy
for insects. If ldr is not available, we use other radar parameters.
Insects are finally screened from liquid layers and melting layer - and
above melting layer.
Args:
obs: The :class:`ClassData` instance.
melting_layer: 2D array denoting melting layer.
liquid_layers: 2D array denoting liquid layers.
prob_lim: Probability higher than this will lead to positive detection.
Default is 0.8.
Returns:
tuple: 2-element tuple containing
- 2-D boolean flag of insects presence.
- 2-D probability of pixel containing insects.
Notes:
This insect detection method is novel and needs to be validated.
"""
probabilities = _insect_probability(obs)
insect_prob = _screen_insects(*probabilities, melting_layer, liquid_layers, obs)
is_insects = insect_prob > prob_lim
return is_insects, ma.masked_where(insect_prob == 0, insect_prob)
def _insect_probability(obs: ClassData) -> tuple[np.ndarray, np.ndarray]:
prob = _get_probabilities(obs)
prob_from_ldr = _calc_prob_from_ldr(prob)
prob_from_others = _calc_prob_from_all(prob)
prob_from_others = _adjust_for_radar(obs, prob, prob_from_others)
prob_combined = _fill_missing_pixels(prob_from_ldr, prob_from_others)
return prob_combined, prob_from_others
def _get_probabilities(obs: ClassData) -> dict:
smooth_v = _get_smoothed_v(obs)
lwp_interp = droplet.interpolate_lwp(obs)
fun = utils.array_to_probability
return {
"width": fun(obs.width, 1, 0.3, invert=True) if hasattr(obs, "width") else 1,
"z_strong": fun(obs.z, 0, 8, invert=True),
"z_weak": fun(obs.z, -20, 8, invert=True),
"ldr": fun(obs.ldr, -25, 5) if hasattr(obs, "ldr") else None,
"sldr": fun(obs.sldr, -25, 5) if hasattr(obs, "sldr") else None,
"temp_loose": fun(obs.tw, 268, 2),
"temp_strict": fun(obs.tw, 274, 1),
"v": fun(smooth_v, -3.5, 2),
"lwp": utils.transpose(fun(lwp_interp, 0.15, 0.05, invert=True)),
"v_sigma": fun(obs.v_sigma, 0.01, 0.1),
}
def _get_smoothed_v(
obs: ClassData,
sigma: tuple[float, float] = (5, 5),
) -> ma.MaskedArray:
smoothed_v = gaussian_filter(obs.v, sigma)
return ma.masked_where(obs.v.mask, smoothed_v)
def _calc_prob_from_ldr(prob: dict) -> np.ndarray:
"""This is the most reliable proxy for insects."""
if prob["ldr"] is not None:
return prob["ldr"] * prob["temp_loose"]
if (
prob["sldr"] is not None
): # Strong SLDR values are probably insects, weak CAN be but not necessarily
p = prob["sldr"]
p[p < 0.9] = ma.masked
return p * prob["temp_loose"]
return np.zeros(prob["z_strong"].shape)
def _calc_prob_from_all(prob: dict) -> np.ndarray:
"""This can be tried when LDR is not available. To detect insects without LDR
unambiguously is difficult and might result in many false positives and/or false
negatives.
"""
return prob["z_weak"] * prob["temp_strict"] * prob["width"] * prob["v"]
def _adjust_for_radar(
obs: ClassData,
prob: dict,
prob_from_others: np.ndarray,
) -> np.ndarray:
"""Adds radar-specific weighting to insect probabilities."""
if "mira" in obs.radar_type.lower():
prob_from_others *= prob["lwp"]
return prob_from_others
def _fill_missing_pixels(
prob_from_ldr: np.ndarray,
prob_from_others: np.ndarray,
) -> np.ndarray:
prob_combined = np.copy(prob_from_ldr)
no_ldr = np.where(prob_from_ldr == 0)
prob_combined[no_ldr] = prob_from_others[no_ldr]
return prob_combined
def _screen_insects(
insect_prob,
insect_prob_no_ldr,
melting_layer,
liquid_layers,
obs,
) -> np.ndarray:
def _screen_liquid_layers() -> None:
prob[liquid_layers == 1] = 0
def _screen_above_melting() -> None:
above_melting = utils.ffill(melting_layer)
prob[above_melting == 1] = 0
def _screen_above_liquid() -> None:
above_liquid = utils.ffill(liquid_layers)
prob[(above_liquid == 1) & (insect_prob_no_ldr > 0)] = 0
def _screen_rainy_profiles() -> None:
rain_smoothed = _smooth_rain(obs.time, obs.is_rain)
prob[rain_smoothed == 1, :] = 0
prob = np.copy(insect_prob)
_screen_liquid_layers()
_screen_above_melting()
_screen_above_liquid()
_screen_rainy_profiles()
return prob
def _smooth_rain(time: np.ndarray, is_rain: np.ndarray) -> np.ndarray:
is_rain_smoothed = np.copy(is_rain)
time_buffer = 5 # minutes
n_profiles = len(is_rain)
n_steps = utils.n_elements(time, time_buffer, "time")
for rain_idx in np.where(is_rain)[0]:
idx_start = max(0, rain_idx - n_steps)
idx_end = min(rain_idx + n_steps, n_profiles)
is_rain_smoothed[idx_start : idx_end + 1] = True
return is_rain_smoothed