Skip to main content

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::field::{ProfileField, WindField};
10use crate::sounding::Sounding;
11
12/// Altitude-layer boundaries, metres above sea level (the sounding's datum,
13/// not above ground): 0–1 km, 1–3 km, 3–6 km, 6–10 km, rest.
14const LAYER_TOPS_M: [f64; 4] = [1_000.0, 3_000.0, 6_000.0, 10_000.0];
15
16/// Mean absolute speed error within one altitude layer.
17#[derive(Debug, Clone, Copy, PartialEq)]
18pub struct LayerError {
19    /// Top of the layer, metres above sea level (`f64::INFINITY` for the
20    /// open layer).
21    pub top_m: f64,
22    /// Mean absolute speed error, m/s.
23    pub mean_abs_speed_error_mps: f64,
24    /// Sounding levels in the layer.
25    pub levels: usize,
26}
27
28/// The comparison result.
29#[derive(Debug, Clone, PartialEq)]
30pub struct ValidationReport {
31    /// Sounding levels compared.
32    pub levels: usize,
33    /// RMS of (profile speed − measured speed), m/s.
34    pub speed_rmse_mps: f64,
35    /// RMS of the circular direction difference, degrees.
36    pub direction_rmse_deg: f64,
37    /// Per-layer mean absolute speed errors.
38    pub layers: Vec<LayerError>,
39}
40
41/// Sample `profile` at every sounding level and compute the metrics. The
42/// profile is read at the sounding's altitudes, above sea level, so its
43/// levels must be in that datum too.
44pub fn validate(profile: &ProfileField, sounding: &Sounding) -> ValidationReport {
45    let mut speed_sq = 0.0;
46    let mut direction_sq = 0.0;
47    let mut layer_abs = vec![0.0; LAYER_TOPS_M.len() + 1];
48    let mut layer_n = vec![0usize; LAYER_TOPS_M.len() + 1];
49    for level in &sounding.levels {
50        let wind = profile.wind(0.0, [0.0, level.altitude_m, 0.0]);
51        let speed_err = wind.speed_mps() - level.speed_mps;
52        speed_sq += speed_err * speed_err;
53        let mut dir_err = (wind.direction_from_deg() - level.direction_from_deg) % 360.0;
54        if dir_err > 180.0 {
55            dir_err -= 360.0;
56        } else if dir_err < -180.0 {
57            dir_err += 360.0;
58        }
59        direction_sq += dir_err * dir_err;
60        let layer = LAYER_TOPS_M
61            .iter()
62            .position(|top| level.altitude_m < *top)
63            .unwrap_or(LAYER_TOPS_M.len());
64        layer_abs[layer] += speed_err.abs();
65        layer_n[layer] += 1;
66    }
67    let n = sounding.levels.len().max(1) as f64;
68    ValidationReport {
69        levels: sounding.levels.len(),
70        speed_rmse_mps: (speed_sq / n).sqrt(),
71        direction_rmse_deg: (direction_sq / n).sqrt(),
72        layers: layer_abs
73            .into_iter()
74            .zip(layer_n)
75            .enumerate()
76            .filter(|(_, (_, count))| *count > 0)
77            .map(|(i, (abs, count))| LayerError {
78                top_m: LAYER_TOPS_M.get(i).copied().unwrap_or(f64::INFINITY),
79                mean_abs_speed_error_mps: abs / count as f64,
80                levels: count,
81            })
82            .collect(),
83    }
84}