sounding_analysis/
wind.rs

1use crate::{
2    error::{AnalysisError, Result},
3    layers::{self, Layer},
4    levels::height_level,
5    sounding::Sounding,
6};
7use itertools::{izip, Itertools};
8use metfor::{IntHelicityM2pS2, Knots, Meters, MetersPSec, Quantity, WindSpdDir, WindUV};
9use std::iter::once;
10
11/// Calculate the mean wind in a layer.
12///
13/// This IS the pressure weighted mean.
14///
15pub fn mean_wind(layer: &Layer, snd: &Sounding) -> Result<WindUV<MetersPSec>> {
16    let pressure = snd.pressure_profile();
17    let wind = snd.wind_profile();
18
19    let top_pres = layer.top.pressure.ok_or(AnalysisError::MissingValue)?;
20    let bottom_pres = layer.bottom.pressure.ok_or(AnalysisError::MissingValue)?;
21    let top_wind = layer.top.wind;
22    let bottom_wind = layer.bottom.wind;
23
24    let intermediate_layers = izip!(pressure, wind)
25        .filter_map(|(pres, wind)| pres.into_option().map(|p| (p, wind)))
26        // Skip values below the layer
27        .skip_while(|&(p, _)| p > bottom_pres)
28        // Only take values below the top of the layer
29        .take_while(|&(p, _)| p > top_pres);
30
31    let (sum_u, sum_v, sum_dp) =
32        // Start at the bottom of the layer
33        once((bottom_pres, &bottom_wind))
34        // Add in any intermediate layers
35        .chain(intermediate_layers)
36        // Finish with the top layer
37        .chain(once((top_pres, &top_wind)))
38        // Filter out missing values
39        .filter_map(|(pres, wind)| wind.map(|w| (pres, w)))
40        // Get the wind u-v components in m/s
41        .map(|(pres, wind)| {
42            let WindUV { u, v } = WindUV::<MetersPSec>::from(wind);
43            (pres, u, v)
44        })
45        // Make windows to see two points at a time for trapezoid rule integration
46        .tuple_windows::<(_, _)>()
47        // Integration with the trapezoid rule to find the mean value
48        .fold(
49            (
50                0.0, // summed u component so far
51                0.0, // summed v component so far
52                0.0, // the total sum of the pressure weights
53            ),
54            |acc, ((p0, u0, v0), (p1, u1, v1))| {
55                let (mut sum_u, mut sum_v, mut acc_dp) = acc;
56
57                let dp = (p0 - p1).unpack();
58
59                sum_u += 0.5 * (u0.unpack() + u1.unpack()) * dp;
60                sum_v += 0.5 * (v0.unpack() + v1.unpack()) * dp;
61                acc_dp += dp;
62
63                (sum_u, sum_v, acc_dp)
64            },
65        );
66
67    if sum_dp == 0.0 {
68        // nothing was done, 1 or zero points in the layer
69        return Err(AnalysisError::NotEnoughData);
70    }
71
72    let avg_u = MetersPSec(sum_u / sum_dp);
73    let avg_v = MetersPSec(sum_v / sum_dp);
74
75    Ok(WindUV { u: avg_u, v: avg_v })
76}
77
78/// Storm relative helicity.
79pub fn sr_helicity<W>(
80    layer: &Layer,
81    storm_motion_uv_ms: W,
82    snd: &Sounding,
83) -> Result<IntHelicityM2pS2>
84where
85    WindUV<MetersPSec>: From<W>,
86{
87    let height = snd.height_profile();
88    let wind = snd.wind_profile();
89    let storm_motion_uv_ms = WindUV::<MetersPSec>::from(storm_motion_uv_ms);
90
91    let bottom = layer.bottom.height.ok_or(AnalysisError::MissingValue)?;
92    let top = layer.top.height.ok_or(AnalysisError::MissingValue)?;
93
94    izip!(height, wind)
95        // Filter out levels with missing values
96        .filter(|(h, w)| h.is_some() && w.is_some())
97        // Unwrap from the `Optioned` type
98        .map(|(h, w)| (h.unpack(), w.unpack()))
99        // Convert the wind and unpack it, subtract the storm motion.
100        .map(|(h, w)| {
101            let WindUV { u, v }: WindUV<MetersPSec> = From::<WindSpdDir<Knots>>::from(w);
102            (h, (u - storm_motion_uv_ms.u), (v - storm_motion_uv_ms.v))
103        })
104        // Make windows so I can see three at a time
105        .tuple_windows::<(_, _, _)>()
106        // Skip levels where the middle of the window is below the bottom of the layer.
107        .skip_while(|(_, (h, _, _), _)| *h < bottom)
108        // Take until the middle of the window is at the top of the layer.
109        .take_while(|(_, (h, _, _), _)| *h <= top)
110        // Add in the derivative information
111        .map(|((h0, u0, v0), (h1, u1, v1), (h2, u2, v2))| {
112            let dz = (h2 - h0).unpack();
113            let du = (u2 - u0).unpack() / dz;
114            let dv = (v2 - v0).unpack() / dz;
115            (h1, u1, v1, du, dv)
116        })
117        // Calculate the helicity (not, this is not the integrated helicity yet)
118        .map(|(z, u, v, du, dv)| (z, u * dv - v * du))
119        // Make a window so I can see two at a time for integrating with the trapezoid rule
120        .tuple_windows::<(_, _)>()
121        // Integrate with the trapezoidal rule
122        .fold(Err(AnalysisError::NotEnoughData), |acc, (lvl0, lvl1)| {
123            let mut integrated_helicity: f64 = acc.unwrap_or(0.0);
124            let (z0, h0) = lvl0;
125            let (z1, h1) = lvl1;
126
127            let h = (h0 + h1).unpack() * (z1 - z0).unpack();
128            integrated_helicity += h;
129
130            Ok(integrated_helicity)
131        })
132        // Multiply by constant factors, -1 for the helicity formula, 2 for trapezoid rule,
133        // and wrap in integrated helicity type
134        .map(|integrated_helicity| IntHelicityM2pS2(-integrated_helicity / 2.0))
135}
136
137/// Calculate the super cell storm motion using the "id" method.
138///
139/// Returns the storm motions in m/s of the right and left mover cells.
140/// (right mover, left mover)
141pub fn bunkers_storm_motion(snd: &Sounding) -> Result<(WindUV<MetersPSec>, WindUV<MetersPSec>)> {
142    let layer = &layers::layer_agl(snd, Meters(6000.0))?;
143
144    let WindUV {
145        u: mean_u,
146        v: mean_v,
147    } = mean_wind(layer, snd)?;
148
149    let WindUV {
150        u: shear_u,
151        v: shear_v,
152    } = bulk_shear_half_km(layer, snd)?;
153
154    const D: f64 = 7.5; // m/s
155
156    let scale = D / shear_u.unpack().hypot(shear_v.unpack());
157    let (delta_u, delta_v) = (shear_v * scale, -shear_u * scale);
158
159    Ok((
160        WindUV {
161            u: mean_u + delta_u,
162            v: mean_v + delta_v,
163        },
164        WindUV {
165            u: mean_u - delta_u,
166            v: mean_v - delta_v,
167        },
168    ))
169}
170
171/// Calculate the bulk shear of a layer using winds averaged over the bottom and top half km.
172///
173/// When using the id method for storm motion vectors, the bulk shear was calculated with top and
174/// bottom wind vectors that were averaged over top/bottom half km of the layer.
175///
176/// Returns `(shear_u_ms, shear_v_ms)`
177pub(crate) fn bulk_shear_half_km(layer: &Layer, snd: &Sounding) -> Result<WindUV<MetersPSec>> {
178    let bottom = layer
179        .bottom
180        .height
181        .into_option()
182        .ok_or(AnalysisError::MissingValue)?;
183    let top = layer
184        .top
185        .height
186        .into_option()
187        .ok_or(AnalysisError::MissingValue)?;
188
189    // abort if not at least 250 meters of non-overlapping area.
190    if top - bottom < Meters(750.0) {
191        return Err(AnalysisError::NotEnoughData);
192    }
193
194    let top_bottom_layer = height_level(bottom + Meters(500.0), snd)?;
195    let bottom_layer = &Layer {
196        top: top_bottom_layer,
197        bottom: layer.bottom,
198    };
199
200    let WindUV {
201        u: bottom_u,
202        v: bottom_v,
203    } = mean_wind(bottom_layer, snd)?;
204
205    let bottom_top_layer = height_level(top - Meters(500.0), snd)?;
206    let top_layer = &Layer {
207        top: layer.top,
208        bottom: bottom_top_layer,
209    };
210
211    let WindUV { u: top_u, v: top_v } = mean_wind(top_layer, snd)?;
212
213    Ok(WindUV {
214        u: top_u - bottom_u,
215        v: top_v - bottom_v,
216    })
217}
sounding_analysis/wind.rs

sounding_analysis/
wind.rs
