winds_engine/
validate.rs

1//! Compare a wind profile against a measured sounding.
2//!
3//! Metrics are sampled at the sounding's own levels: speed RMSE, direction
4//! RMSE (circular), and mean absolute speed error per altitude layer
5//! (1, 3, 6 and 10 km boundaries, then the rest). Layers bin on the
6//! sounding's altitudes, which are above sea level, not above ground: for a
7//! station at 287 m the first layer spans only 713 m above ground.
8
9use crate::sounding::Sounding;
10use crate::ProfileField;
11use nalgebra::Vector3;
12
13/// Horizontal speed of a wind vector, m/s.
14fn horizontal_speed(w: Vector3<f64>) -> f64 {
15    w.x.hypot(w.z)
16}
17
18/// Wind-FROM direction of a wind vector, degrees in [0, 360) (0° = headwind,
19/// 90° = from the right).
20fn direction_from_deg(w: Vector3<f64>) -> f64 {
21    (w.z.atan2(w.x).to_degrees() + 180.0).rem_euclid(360.0)
22}
23
24/// Altitude-layer boundaries, metres above sea level (the sounding's datum,
25/// not above ground): 0–1 km, 1–3 km, 3–6 km, 6–10 km, rest.
26const LAYER_TOPS_M: [f64; 4] = [1_000.0, 3_000.0, 6_000.0, 10_000.0];
27
28/// Mean absolute speed error within one altitude layer.
29#[derive(Debug, Clone, Copy, PartialEq)]
30pub struct LayerError {
31    /// Top of the layer, metres above sea level (`f64::INFINITY` for the
32    /// open layer).
33    pub top_m: f64,
34    /// Mean absolute speed error, m/s.
35    pub mean_abs_speed_error_mps: f64,
36    /// Sounding levels in the layer.
37    pub levels: usize,
38}
39
40/// The comparison result.
41#[derive(Debug, Clone, PartialEq)]
42pub struct ValidationReport {
43    /// Sounding levels compared.
44    pub levels: usize,
45    /// RMS of (profile speed − measured speed), m/s.
46    pub speed_rmse_mps: f64,
47    /// RMS of the circular direction difference, degrees.
48    pub direction_rmse_deg: f64,
49    /// Per-layer mean absolute speed errors.
50    pub layers: Vec<LayerError>,
51}
52
53/// Sample `profile` at every sounding level and compute the metrics. The
54/// profile is read at the sounding's altitudes, above sea level, so its
55/// levels must be in that datum too.
56pub fn validate(profile: &ProfileField, sounding: &Sounding) -> ValidationReport {
57    let mut speed_sq = 0.0;
58    let mut direction_sq = 0.0;
59    let mut layer_abs = vec![0.0; LAYER_TOPS_M.len() + 1];
60    let mut layer_n = vec![0usize; LAYER_TOPS_M.len() + 1];
61    for level in &sounding.levels {
62        let wind = profile.wind_at(level.altitude_m);
63        let speed_err = horizontal_speed(wind) - level.speed_mps;
64        speed_sq += speed_err * speed_err;
65        let mut dir_err = (direction_from_deg(wind) - level.direction_from_deg) % 360.0;
66        if dir_err > 180.0 {
67            dir_err -= 360.0;
68        } else if dir_err < -180.0 {
69            dir_err += 360.0;
70        }
71        direction_sq += dir_err * dir_err;
72        let layer = LAYER_TOPS_M
73            .iter()
74            .position(|top| level.altitude_m < *top)
75            .unwrap_or(LAYER_TOPS_M.len());
76        layer_abs[layer] += speed_err.abs();
77        layer_n[layer] += 1;
78    }
79    let n = sounding.levels.len().max(1) as f64;
80    ValidationReport {
81        levels: sounding.levels.len(),
82        speed_rmse_mps: (speed_sq / n).sqrt(),
83        direction_rmse_deg: (direction_sq / n).sqrt(),
84        layers: layer_abs
85            .into_iter()
86            .zip(layer_n)
87            .enumerate()
88            .filter(|(_, (_, count))| *count > 0)
89            .map(|(i, (abs, count))| LayerError {
90                top_m: LAYER_TOPS_M.get(i).copied().unwrap_or(f64::INFINITY),
91                mean_abs_speed_error_mps: abs / count as f64,
92                levels: count,
93            })
94            .collect(),
95    }
96}
winds_engine/validate.rs

winds_engine/
validate.rs
