sounding_analysis/
levels.rs

1//! This module finds significant levels such as the freezing level and wet bulb zero level. It also
2//! has functions for finding critical values at a single level, such as the maximum wet bulb
3//! temperature aloft.  It does not include functions for finding levels related to parcel analysis
4//! and convection, those are found in the `parcel` module.
5use crate::{
6    error::{
7        AnalysisError::{InvalidInput, MissingProfile},
8        {AnalysisError, Result},
9    },
10    interpolation::{linear_interp, linear_interpolate_sounding},
11    layers::Layer,
12    sounding::{DataRow, Sounding},
13};
14use itertools::{izip, Itertools};
15use metfor::{Celsius, HectoPascal, Meters, FREEZING};
16use optional::Optioned;
17
18/// A level in the atmosphere is described by a `DataRow` from a sounding.
19pub type Level = DataRow;
20
21/// A list of levels.
22pub type Levels = Vec<Level>;
23
24/// Find the freezing/melting levels below 500 hPa.
25pub fn freezing_levels(snd: &Sounding) -> Result<Levels> {
26    find_temperature_levels(
27        FREEZING,
28        snd.pressure_profile(),
29        snd.temperature_profile(),
30        snd,
31    )
32}
33
34/// Find the wet bulb zero levels
35pub fn wet_bulb_zero_levels(snd: &Sounding) -> Result<Levels> {
36    find_temperature_levels(
37        FREEZING,
38        snd.pressure_profile(),
39        snd.wet_bulb_profile(),
40        snd,
41    )
42}
43
44fn find_temperature_levels(
45    target_t: Celsius,
46    p_profile: &[Optioned<HectoPascal>],
47    t_profile: &[Optioned<Celsius>],
48    snd: &Sounding,
49) -> Result<Levels> {
50    const TOP_PRESSURE: HectoPascal = HectoPascal(500.0); // don't look above here
51
52    if t_profile.is_empty() || p_profile.is_empty() {
53        return Err(AnalysisError::MissingProfile);
54    }
55
56    #[inline]
57    fn crosses_target(t0: Celsius, t1: Celsius, target_t: Celsius) -> bool {
58        (t0 <= target_t && t1 >= target_t) || (t0 >= target_t && t1 <= target_t)
59    }
60
61    izip!(p_profile, t_profile)
62        // Remove levels with missing data
63        .filter(|(p, t)| p.is_some() && t.is_some())
64        // Unwrap from the Optioned type
65        .map(|(p, t)| (p.unpack(), t.unpack()))
66        // Don't bother looking above a certain level
67        .take_while(|&(p, _)| p >= TOP_PRESSURE)
68        // Look at them in pairs to find crossing the threshold
69        .tuple_windows::<(_, _)>()
70        // Filter out any pair that is not crossing the desired temperature level
71        .filter(|&((_p0, t0), (_p1, t1))| crosses_target(t0, t1, target_t))
72        // Interpolate crossings to get pressures values at the levels we want.
73        .map(|((p0, t0), (p1, t1))| linear_interp(target_t, t0, t1, p0, p1))
74        // Perform the interpolation on the full sounding to get the desired DataRow
75        .map(|p_level| linear_interpolate_sounding(snd, p_level))
76        .collect()
77}
78
79/// Maximum wet bulb temperature aloft.
80pub fn max_wet_bulb_in_profile(snd: &Sounding) -> Result<Level> {
81    max_t_aloft(snd, snd.wet_bulb_profile())
82}
83
84/// Maximum temperature aloft.
85pub fn max_temperature_in_profile(snd: &Sounding) -> Result<Level> {
86    max_t_aloft(snd, snd.temperature_profile())
87}
88
89// Only searches up to 500 hPa
90fn max_t_aloft(snd: &Sounding, t_profile: &[Optioned<Celsius>]) -> Result<Level> {
91    const TOP_PRESSURE: HectoPascal = HectoPascal(500.0); // don't look above here.
92
93    let p_profile = snd.pressure_profile();
94
95    if t_profile.is_empty() || p_profile.is_empty() {
96        return Err(AnalysisError::MissingProfile);
97    }
98
99    let sentinel = Celsius::from(metfor::Kelvin(0.0));
100
101    let (idx, mx_t) = izip!(0usize.., p_profile, t_profile)
102        // Filter out levels with missing values
103        .filter(|(_, p, t)| p.is_some() && t.is_some())
104        // unwrap the Optioned type
105        .map(|(i, p, t)| (i, p.unpack(), t.unpack()))
106        // Only look up to a certain level
107        .take_while(|&(_, p, _)| p >= TOP_PRESSURE)
108        // fold to get the maximum value
109        .fold((0, sentinel), |acc, (i, _, t)| {
110            let (_, mx_t) = acc;
111            if t > mx_t {
112                (i, t)
113            } else {
114                // Propagate most recent result through
115                acc
116            }
117        });
118
119    if mx_t == sentinel {
120        Err(AnalysisError::NotEnoughData)
121    } else {
122        snd.data_row(idx).ok_or(InvalidInput)
123    }
124}
125
126/// Maximum temperature in a layer.
127pub fn max_temperature_in_layer(snd: &Sounding, lyr: &Layer) -> Result<Level> {
128    max_t_in_layer(snd, snd.temperature_profile(), lyr)
129}
130
131/// Maximum wet bulb temperature in a layer.
132pub fn max_wet_bulb_in_layer(snd: &Sounding, lyr: &Layer) -> Result<Level> {
133    max_t_in_layer(snd, snd.wet_bulb_profile(), lyr)
134}
135
136#[inline]
137fn max_t_in_layer(snd: &Sounding, t_profile: &[Optioned<Celsius>], lyr: &Layer) -> Result<Level> {
138    let (bottom_p, top_p) = if lyr.bottom.pressure.is_some() && lyr.top.pressure.is_some() {
139        (lyr.bottom.pressure.unpack(), lyr.top.pressure.unpack())
140    } else {
141        return Err(AnalysisError::InvalidInput);
142    };
143
144    let p_profile = snd.pressure_profile();
145
146    if t_profile.is_empty() || p_profile.is_empty() {
147        return Err(AnalysisError::MissingProfile);
148    }
149
150    izip!(0usize.., p_profile, t_profile)
151        // Remove levels with missing values
152        .filter(|(_, p, t)| p.is_some() && t.is_some())
153        // Unwrap the Optioned type
154        .map(|(i, p, t)| (i, p.unpack(), t.unpack()))
155        // Skip values below the bottom of the layer
156        .skip_while(|(_, p, _)| *p > bottom_p)
157        // Only look up to the top of the layer
158        .take_while(|(_, p, _)| *p >= top_p)
159        // fold to find the max value
160        .fold(
161            Err(AnalysisError::NotEnoughData),
162            |acc: Result<_>, (i, _, t)| {
163                if let Ok((_, mx_t)) = acc {
164                    if t > mx_t {
165                        Ok((i, t))
166                    } else {
167                        // Propagate most recent result through
168                        acc
169                    }
170                } else {
171                    // We're just starting, so populate the result
172                    Ok((i, t))
173                }
174            },
175        )
176        // Retrive the row
177        .and_then(|(idx, _)| snd.data_row(idx).ok_or(InvalidInput))
178}
179
180/// Find a level at a specific geopotential height.
181pub(crate) fn height_level(tgt_height: Meters, snd: &Sounding) -> Result<Level> {
182    let h_profile = snd.height_profile();
183    let p_profile = snd.pressure_profile();
184
185    if h_profile.is_empty() || p_profile.is_empty() {
186        return Err(MissingProfile);
187    }
188
189    // Assumes h0 < h1, ie iterate from bottom to top of sounding
190    #[inline]
191    fn is_cross_over(h0: Meters, h1: Meters, tgt_height: Meters) -> bool {
192        debug_assert!(h0 <= h1);
193        h0 <= tgt_height && h1 >= tgt_height
194    }
195
196    izip!(p_profile, h_profile)
197        // Filter out levels with missing data
198        .filter(|(p, h)| p.is_some() && h.is_some())
199        // Unpack from the Optioned type
200        .map(|(p, h)| (p.unpack(), h.unpack()))
201        // look at levels in pairs to find when the data crosses the desired height
202        .tuple_windows::<(_, _)>()
203        // Filter out all pairs that don't have a crossover
204        .find(|&((_p0, h0), (_p1, h1))| is_cross_over(h0, h1, tgt_height))
205        // Map the bracket into the pressure at the target height via interpolation
206        .map(|((p0, h0), (p1, h1))| linear_interp(tgt_height, h0, h1, p0, p1))
207        // map the option into a result
208        .ok_or(AnalysisError::InvalidInput)
209        // Interpolate the result
210        .and_then(|level_p| linear_interpolate_sounding(snd, level_p))
211}
sounding_analysis/levels.rs

sounding_analysis/
levels.rs
