use crate::{
error::{
AnalysisError::{InvalidInput, MissingProfile},
{AnalysisError, Result},
},
interpolation::{linear_interp, linear_interpolate_sounding},
layers::Layer,
sounding::{DataRow, Sounding},
};
use itertools::{izip, Itertools};
use metfor::{Celsius, HectoPascal, Meters, FREEZING};
use optional::Optioned;
pub type Level = DataRow;
pub type Levels = Vec<Level>;
pub fn freezing_levels(snd: &Sounding) -> Result<Levels> {
find_temperature_levels(
FREEZING,
snd.pressure_profile(),
snd.temperature_profile(),
snd,
)
}
pub fn wet_bulb_zero_levels(snd: &Sounding) -> Result<Levels> {
find_temperature_levels(
FREEZING,
snd.pressure_profile(),
snd.wet_bulb_profile(),
snd,
)
}
fn find_temperature_levels(
target_t: Celsius,
p_profile: &[Optioned<HectoPascal>],
t_profile: &[Optioned<Celsius>],
snd: &Sounding,
) -> Result<Levels> {
const TOP_PRESSURE: HectoPascal = HectoPascal(500.0);
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
#[inline]
fn crosses_target(t0: Celsius, t1: Celsius, target_t: Celsius) -> bool {
(t0 <= target_t && t1 >= target_t) || (t0 >= target_t && t1 <= target_t)
};
izip!(p_profile, t_profile)
.filter(|(p, t)| p.is_some() && t.is_some())
.map(|(p, t)| (p.unpack(), t.unpack()))
.take_while(|&(p, _)| p >= TOP_PRESSURE)
.tuple_windows::<(_, _)>()
.filter(|&((_p0, t0), (_p1, t1))| crosses_target(t0, t1, target_t))
.map(|((p0, t0), (p1, t1))| linear_interp(target_t, t0, t1, p0, p1))
.map(|p_level| linear_interpolate_sounding(snd, p_level))
.collect()
}
pub fn max_wet_bulb_in_profile(snd: &Sounding) -> Result<Level> {
max_t_aloft(snd, snd.wet_bulb_profile())
}
pub fn max_temperature_in_profile(snd: &Sounding) -> Result<Level> {
max_t_aloft(snd, snd.temperature_profile())
}
fn max_t_aloft(snd: &Sounding, t_profile: &[Optioned<Celsius>]) -> Result<Level> {
const TOP_PRESSURE: HectoPascal = HectoPascal(500.0);
let p_profile = snd.pressure_profile();
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
izip!(0usize.., p_profile, t_profile)
.filter(|(_, p, t)| p.is_some() && t.is_some())
.map(|(i, p, t)| (i, p.unpack(), t.unpack()))
.take_while(|&(_, p, _)| p >= TOP_PRESSURE)
.fold(
Err(AnalysisError::NotEnoughData),
|acc: Result<_>, (i, _, t)| {
if let Ok((_, mx_t)) = acc {
if t > mx_t {
Ok((i, t))
} else {
acc
}
} else {
Ok((i, t))
}
},
)
.and_then(|(idx, _)| snd.data_row(idx).ok_or(InvalidInput))
}
pub fn max_temperature_in_layer(snd: &Sounding, lyr: &Layer) -> Result<Level> {
max_t_in_layer(snd, snd.temperature_profile(), lyr)
}
pub fn max_wet_bulb_in_layer(snd: &Sounding, lyr: &Layer) -> Result<Level> {
max_t_in_layer(snd, snd.wet_bulb_profile(), lyr)
}
#[inline]
fn max_t_in_layer(snd: &Sounding, t_profile: &[Optioned<Celsius>], lyr: &Layer) -> Result<Level> {
let (bottom_p, top_p) = if lyr.bottom.pressure.is_some() && lyr.top.pressure.is_some() {
(lyr.bottom.pressure.unpack(), lyr.top.pressure.unpack())
} else {
return Err(AnalysisError::InvalidInput);
};
let p_profile = snd.pressure_profile();
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
izip!(0usize.., p_profile, t_profile)
.filter(|(_, p, t)| p.is_some() && t.is_some())
.map(|(i, p, t)| (i, p.unpack(), t.unpack()))
.skip_while(|(_, p, _)| *p > bottom_p)
.take_while(|(_, p, _)| *p >= top_p)
.fold(
Err(AnalysisError::NotEnoughData),
|acc: Result<_>, (i, _, t)| {
if let Ok((_, mx_t)) = acc {
if t > mx_t {
Ok((i, t))
} else {
acc
}
} else {
Ok((i, t))
}
},
)
.and_then(|(idx, _)| snd.data_row(idx).ok_or(InvalidInput))
}
pub(crate) fn height_level(tgt_height: Meters, snd: &Sounding) -> Result<Level> {
let h_profile = snd.height_profile();
let p_profile = snd.pressure_profile();
if h_profile.is_empty() || p_profile.is_empty() {
return Err(MissingProfile);
}
#[inline]
fn is_cross_over(h0: Meters, h1: Meters, tgt_height: Meters) -> bool {
debug_assert!(h0 <= h1);
h0 <= tgt_height && h1 >= tgt_height
}
izip!(p_profile, h_profile)
.filter(|(p, h)| p.is_some() && h.is_some())
.map(|(p, h)| (p.unpack(), h.unpack()))
.tuple_windows::<(_, _)>()
.filter(|&((_p0, h0), (_p1, h1))| is_cross_over(h0, h1, tgt_height))
.nth(0)
.map(|((p0, h0), (p1, h1))| linear_interp(tgt_height, h0, h1, p0, p1))
.ok_or(AnalysisError::InvalidInput)
.and_then(|level_p| linear_interpolate_sounding(snd, level_p))
}