sounding_analysis/layers/
inversions.rs

1use super::{Layer, Layers};
2use crate::{
3    error::{AnalysisError, Result},
4    sounding::Sounding,
5};
6use itertools::{izip, Itertools};
7use metfor::HectoPascal;
8
9/// Get all inversion layers up to a specified pressure.
10pub fn inversions(snd: &Sounding, top_p: HectoPascal) -> Result<Layers> {
11    let t_profile = snd.temperature_profile();
12    let p_profile = snd.pressure_profile();
13
14    if t_profile.is_empty() || p_profile.is_empty() {
15        return Err(AnalysisError::MissingProfile);
16    }
17
18    izip!(0usize.., p_profile, t_profile)
19        // Filter out levels with missing data
20        .filter(|&(_, p, t)| p.is_some() && t.is_some())
21        // Unpack from the `Optioned` type
22        .map(|(i, p, t)| (i, p.unpack(), t.unpack()))
23        // Stop once we reach the top level pressure
24        .take_while(|&(_, p, _)| p >= top_p)
25        // Now all we need is the row index and the temperature
26        .map(|(i, _, t)| (i, t))
27        // Inspect the levels in pairs
28        .tuple_windows::<(_, _)>()
29        // Map into inversion types
30        .map(|((i0, t0), (i1, t1))| {
31            if t1 > t0 {
32                LayerType::Inversion(i0, i1)
33            } else {
34                LayerType::NonInversion
35            }
36        })
37        // Combine adjacent inversion layers into a single, larger inversion
38        .scan(
39            (None, None),
40            |layers: &mut (Option<usize>, Option<usize>), layer: LayerType| {
41                let bottom: &mut Option<usize> = &mut layers.0;
42                let top: &mut Option<usize> = &mut layers.1;
43
44                match layer {
45                    LayerType::Inversion(i0, i1) => {
46                        if bottom.is_none() {
47                            *bottom = Some(i0);
48                        }
49                        *top = Some(i1);
50                        Some(None)
51                    }
52                    LayerType::NonInversion => {
53                        if let (Some(i0), Some(i1)) = (bottom.take(), top.take()) {
54                            Some(Some((i0, i1)))
55                        } else {
56                            Some(None)
57                        }
58                    }
59                }
60            },
61        )
62        // Remove iterations that did not capture an inverison
63        .flatten()
64        // Map indexes into layers
65        .map(|(i0, i1)| {
66            snd.data_row(i0)
67                .and_then(|bottom| snd.data_row(i1).map(|top| Layer { bottom, top }))
68        })
69        // Transform from `Option` to `Result`
70        .map(|opt| opt.ok_or(AnalysisError::InvalidInput))
71        .collect()
72}
73
74/// Get a surface based inversion.
75pub fn sfc_based_inversion(snd: &Sounding) -> Result<Option<Layer>> {
76    let t_profile = snd.temperature_profile();
77    let p_profile = snd.pressure_profile();
78    let h_profile = snd.height_profile();
79
80    if t_profile.is_empty() || p_profile.is_empty() || h_profile.is_empty() {
81        return Err(AnalysisError::MissingProfile);
82    }
83
84    let sfc_layer = izip!(0usize.., p_profile, h_profile, t_profile)
85        // Remove levels with missing data
86        .filter(|(_, p, h, t)| p.is_some() && h.is_some() && t.is_some())
87        // Unpack from `Optioned` type and discard the height, we no longer need it.
88        .map(|(i, p, _, t)| (i, p.unpack(), t.unpack()))
89        // Only look up to about 700 hPa
90        .take_while(|(_, p, _)| *p > HectoPascal(690.0))
91        // Pair them up to look at them as a layer
92        .tuple_windows::<(_, _)>()
93        // Map into inversion types
94        .map(|((i0, _, t0), (i1, _, t1))| {
95            if t1 > t0 {
96                LayerType::Inversion(i0, i1)
97            } else {
98                LayerType::NonInversion
99            }
100        })
101        // Only take inversion layers to be combined later, if the first one isn't an inversion,
102        // then this isn't a surface based inversion!
103        .take_while(|lyr| {
104            std::mem::discriminant(lyr) == std::mem::discriminant(&LayerType::Inversion(0, 0))
105        })
106        // Unwrap from the inversion type
107        .map(|lyr| match lyr {
108            LayerType::Inversion(i0, i1) => (i0, i1),
109            LayerType::NonInversion => unreachable!(),
110        })
111        // Now combine the layers into one
112        .fold(None, |combined_layers, (i0, i1)| match combined_layers {
113            None => Some((i0, i1)),
114            Some((bottom, _top)) => Some((bottom, i1)),
115        }) // Option<(usize, usize)
116        .and_then(|(bottom_i, top_i)| {
117            let bottom = snd.data_row(bottom_i)?;
118            let top = snd.data_row(top_i)?;
119            Some(Layer { bottom, top })
120        });
121
122    Ok(sfc_layer)
123}
124
125enum LayerType {
126    Inversion(usize, usize),
127    NonInversion,
128}
sounding_analysis/layers/inversions.rs

sounding_analysis/layers/
inversions.rs
