sounding_analysis/
profile.rs

1//! Create profiles.
2//!
3//! The output will be at the same levels as the sounding and these are suitable to be set as a
4//! profile in the sounding. For example, calculating a wet bulb or relative humidity profile from a
5//! sounding with temperature and dew point. If one of the profiles required for the analysis in the
6//! sounding is missing, the result cannot be calculated and an empty vector is returned.
7//!
8use crate::sounding::Sounding;
9use itertools::izip;
10use metfor::{
11    self, Celsius, CelsiusDiff, CelsiusPKm, HydrolapsePKm, Kelvin, KelvinPKm, Km, Meters, Quantity,
12    Temperature,
13};
14use optional::{none, some, Optioned};
15use std::ops::Sub;
16
17/// Given a sounding, calculate a profile of wet bulb temperature.
18pub fn wet_bulb(snd: &Sounding) -> Vec<Optioned<Celsius>> {
19    let p_profile = snd.pressure_profile();
20    let t_profile = snd.temperature_profile();
21    let dp_profile = snd.dew_point_profile();
22
23    if p_profile.len().min(t_profile.len()).min(dp_profile.len()) == 0 {
24        return vec![];
25    }
26
27    izip!(p_profile, t_profile, dp_profile)
28        .map(|(p_opt, t_opt, dp_opt)| {
29            p_opt.and_then(|p| {
30                t_opt.and_then(|t| {
31                    dp_opt.and_then(|dp| Optioned::<Celsius>::from(metfor::wet_bulb(t, dp, p)))
32                })
33            })
34        })
35        .collect()
36}
37
38/// Given a sounding, calculate a profile of relative humidity.
39pub fn relative_humidity(snd: &Sounding) -> Vec<Optioned<f64>> {
40    let t_profile = snd.temperature_profile();
41    let dp_profile = snd.dew_point_profile();
42
43    if t_profile.len().min(dp_profile.len()) == 0 {
44        return vec![];
45    }
46
47    izip!(t_profile, dp_profile)
48        .map(|(t_opt, dp_opt)| {
49            t_opt.and_then(|t| dp_opt.and_then(|dp| Optioned::from(metfor::rh(t, dp))))
50        })
51        .collect()
52}
53
54/// Given a sounding, calculate a profile of relative humidity with respect to ice.
55pub fn relative_humidity_ice(snd: &Sounding) -> Vec<Optioned<f64>> {
56    let t_profile = snd.temperature_profile();
57    let dp_profile = snd.dew_point_profile();
58
59    if t_profile.len().min(dp_profile.len()) == 0 {
60        return vec![];
61    }
62
63    izip!(t_profile, dp_profile)
64        .map(|(t_opt, dp_opt)| {
65            t_opt.and_then(|t| {
66                dp_opt.and_then(|dp| {
67                    // Use dew point to get vapor pressure of water - sounding more likely to have
68                    // dew point information.
69                    let vp_water = metfor::vapor_pressure_water(dp);
70                    // Get the saturation vapor pressure relative to ice
71                    let vp_sat_ice = metfor::vapor_pressure_ice(t);
72
73                    let rh = vp_water.and_then(|vpw| vp_sat_ice.map(|vpsat| vpw / vpsat));
74                    Optioned::from(rh)
75                })
76            })
77        })
78        .collect()
79}
80
81/// Given a sounding, calculate a profile of the potential temperature.
82pub fn potential_temperature(snd: &Sounding) -> Vec<Optioned<Kelvin>> {
83    let p_profile = snd.pressure_profile();
84    let t_profile = snd.temperature_profile();
85
86    if p_profile.len().min(t_profile.len()) == 0 {
87        return vec![];
88    }
89
90    izip!(p_profile, t_profile)
91        .map(|(p_opt, t_opt)| {
92            p_opt.and_then(|p| {
93                t_opt.and_then(|t| Optioned::<Kelvin>::from(metfor::potential_temperature(p, t)))
94            })
95        })
96        .collect()
97}
98
99/// Given a sounding, calculate a profile of the equivalent potential temperature.
100pub fn equivalent_potential_temperature(snd: &Sounding) -> Vec<Optioned<Kelvin>> {
101    let p_profile = snd.pressure_profile();
102    let t_profile = snd.temperature_profile();
103    let dp_profile = snd.dew_point_profile();
104
105    if p_profile.len().min(t_profile.len()).min(dp_profile.len()) == 0 {
106        return vec![];
107    }
108
109    izip!(p_profile, t_profile, dp_profile)
110        .map(|(p_opt, t_opt, dp_opt)| {
111            p_opt.and_then(|p| {
112                t_opt.and_then(|t| {
113                    dp_opt.and_then(|dp| {
114                        Optioned::<Kelvin>::from(metfor::equiv_pot_temperature(t, dp, p))
115                    })
116                })
117            })
118        })
119        .collect()
120}
121
122/// Get a profile of the lapse rate between layers in &deg;C / km.
123pub fn temperature_lapse_rate(snd: &Sounding) -> Vec<Optioned<CelsiusPKm>> {
124    let t_profile = snd.temperature_profile().iter().cloned();
125    lapse_rate(snd, t_profile)
126}
127
128/// Get a profile of the average lapse rate from the surface to *, or the level on the y axis.
129pub fn sfc_to_level_temperature_lapse_rate(snd: &Sounding) -> Vec<Optioned<CelsiusPKm>> {
130    let z_profile = snd.height_profile();
131    let t_profile = snd.temperature_profile();
132
133    let (t_sfc, z_sfc): (Celsius, Meters) = if let (Some(t_sfc), Some(z_sfc)) = (
134        snd.sfc_temperature().into(),
135        snd.station_info().elevation().into(),
136    ) {
137        (t_sfc, z_sfc)
138    } else {
139        return vec![];
140    };
141
142    izip!(z_profile, t_profile)
143        .map(|(z_opt, t_opt)| {
144            if let (Some(z), Some(t)) = (z_opt.into_option(), t_opt.into_option()) {
145                if (z - z_sfc).abs() < Meters(std::f64::EPSILON) {
146                    none()
147                } else {
148                    let CelsiusDiff(dt) = t - t_sfc;
149                    let Km(dz) = Km::from(z - z_sfc);
150                    some(CelsiusPKm(dt / dz))
151                }
152            } else {
153                none()
154            }
155        })
156        .collect()
157}
158
159/// Get the lapse rate of equivalent potential temperature in &deg;K / km.
160pub fn theta_e_lapse_rate(snd: &Sounding) -> Vec<Optioned<KelvinPKm>> {
161    let theta_e = snd.theta_e_profile().iter().cloned();
162    lapse_rate(snd, theta_e)
163}
164
165fn lapse_rate<I, T>(snd: &Sounding, v_profile: I) -> Vec<Optioned<CelsiusPKm>>
166where
167    I: Iterator<Item = Optioned<T>>,
168    T: Temperature + Sub<T> + optional::Noned,
169    CelsiusDiff: From<<T as Sub<T>>::Output>,
170{
171    let z_profile = snd.height_profile();
172
173    izip!(z_profile, v_profile)
174        .scan((None, None), |prev_pair, (&z, v)| {
175            let &mut (ref mut prev_z, ref mut prev_v) = prev_pair;
176
177            let z: Option<Meters> = z.into_option();
178            let v: Option<T> = v.into_option();
179
180            let lapse_rate = if let (Some(ref prev_z), Some(ref prev_v), Some(ref z), Some(ref v)) =
181                (*prev_z, *prev_v, z, v)
182            {
183                let CelsiusDiff(dt) = CelsiusDiff::from(*v - *prev_v);
184                let Km(dz) = Km::from(*z - *prev_z);
185                some(CelsiusPKm(dt / dz))
186            } else {
187                none()
188            };
189
190            *prev_z = z;
191            *prev_v = v;
192
193            Some(lapse_rate)
194        })
195        .collect()
196}
197
198/// Get the hydrolapse in (kg/kg)/km
199pub fn hydrolapse(snd: &Sounding) -> Vec<Optioned<HydrolapsePKm>> {
200    let z_profile = snd.height_profile();
201    let dp_profile = snd.dew_point_profile();
202    let p_profile = snd.pressure_profile();
203
204    izip!(p_profile, z_profile, dp_profile)
205        .scan((None, None), |prev_pair, (&p, &z, &dp)| {
206            let &mut (ref mut prev_z, ref mut prev_mw) = prev_pair;
207
208            let p: Option<_> = p.into_option();
209            let z: Option<_> = z.into_option();
210            let dp: Option<_> = dp.into_option();
211
212            let mw = if let (Some(p), Some(dp)) = (p, dp) {
213                metfor::mixing_ratio(dp, p)
214            } else {
215                None
216            };
217
218            let mw_lapse_rate =
219                if let (Some(p_z), Some(p_mw), Some(z), Some(mw)) = (*prev_z, *prev_mw, z, mw) {
220                    let dmw = mw - p_mw;
221                    let Km(dz) = Km::from(z - p_z);
222                    some(HydrolapsePKm(dmw / dz))
223                } else {
224                    none()
225                };
226
227            *prev_z = z;
228            *prev_mw = mw;
229
230            Some(mw_lapse_rate)
231        })
232        .collect()
233}
234
235#[cfg(test)]
236mod test_sounding_profiles {
237    use super::*;
238    use metfor::{Celsius, Meters};
239
240    fn make_test_sounding() -> Sounding {
241        Sounding::new()
242            .with_temperature_profile(vec![
243                some(Celsius(9.8)),
244                some(Celsius(0.0)),
245                some(Celsius(-5.0)),
246            ])
247            .with_height_profile(vec![
248                some(Meters(1000.0)),
249                some(Meters(2000.0)),
250                some(Meters(3000.0)),
251            ])
252    }
253
254    #[test]
255    fn test_temperature_lapse_rate() {
256        let snd = make_test_sounding();
257
258        let lapse_rate = temperature_lapse_rate(&snd);
259        println!("{:#?}", lapse_rate);
260        assert!(lapse_rate.contains(&some(CelsiusPKm(-9.8))));
261        assert!(lapse_rate.contains(&some(CelsiusPKm(-5.0))));
262    }
263}
sounding_analysis/profile.rs

sounding_analysis/
profile.rs
