Graham McIntire
af560be01e
Clippy's `neg_cmp_op_on_partial_ord` lint fails on `-D warnings`, which CI enforces. Using `<=` directly is clearer and equivalent for the f64 inputs we pass (both are caller-supplied target pressures that come from the level-key constants — neither is ever NaN).
526 lines
18 KiB
Rust
526 lines
18 KiB
Rust
//! Per-cell *derived* weather scalars on disk, the cheap-read sibling
|
||
//! of [`profiles_file`][crate::profiles_file].
|
||
//!
|
||
//! 1:1 wire-compatible with `Microwaveprop.Weather.ScalarFile` on the
|
||
//! Elixir side. Chunked into 5°×5° spatial buckets so a state-sized
|
||
//! `/weather` viewport reads only the chunks that overlap.
|
||
//!
|
||
//! ## Layout
|
||
//!
|
||
//! ```text
|
||
//! <scores_dir>/weather_scalars/
|
||
//! <iso>/ # e.g. 2026-04-29T12:00:00Z/
|
||
//! <lat_band>_<lon_band>.mp.gz
|
||
//! ```
|
||
//!
|
||
//! `lat_band = floor(lat / 5)`, `lon_band = floor(lon / 5)`. Each chunk
|
||
//! is gzipped MessagePack: `[ {row}, {row}, … ]`. Map keys are strings
|
||
//! on the wire — the Elixir reader atomizes the whitelist on read, so
|
||
//! every key emitted here must match the Elixir whitelist or it gets
|
||
//! silently dropped by callers that use atom-style access.
|
||
//!
|
||
//! Scalar derivation mirrors `Microwaveprop.Weather.WeatherLayers.derive/1`
|
||
//! plus the surface fields `Microwaveprop.Weather.build_grid_cache_row/4`
|
||
//! adds on top. The Elixir module remains the single source of truth
|
||
//! for any cell that hasn't been Rust-materialized yet (the
|
||
//! `kickoff_async_scalar_materialize` fallback path).
|
||
|
||
use std::collections::HashMap;
|
||
use std::fs::File;
|
||
use std::io::{BufWriter, Write};
|
||
use std::path::{Path, PathBuf};
|
||
|
||
use chrono::{DateTime, Utc};
|
||
use flate2::write::GzEncoder;
|
||
use flate2::Compression;
|
||
use serde::Serialize;
|
||
|
||
use crate::decoder::CellValues;
|
||
use crate::fetcher;
|
||
use crate::sounding_params::{
|
||
ducting_detected, min_refractivity_gradient, sat_vap_pres, surface_refractivity, Level,
|
||
};
|
||
|
||
const CHUNK_STEP: i32 = 5;
|
||
const SUBDIR: &str = "weather_scalars";
|
||
|
||
/// One row in a chunk file. Field names match the Elixir atom-key
|
||
/// whitelist exactly — adding a field here without updating the
|
||
/// Elixir `@atom_keys` set means the new field round-trips as a
|
||
/// string and is invisible to atom-style callers.
|
||
#[derive(Debug, Default, Serialize)]
|
||
pub struct ScalarRow {
|
||
pub lat: f64,
|
||
pub lon: f64,
|
||
pub valid_time: String,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub temperature: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub dewpoint_depression: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub surface_rh: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub surface_pressure_mb: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub surface_refractivity: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub refractivity_gradient: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub bl_height: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub pwat: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub temp_850mb: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub dewpoint_850mb: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub temp_700mb: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub dewpoint_700mb: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub lapse_rate: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub mid_lapse_rate: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub inversion_strength: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub inversion_base_m: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub ducting: Option<bool>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub duct_base_m: Option<f64>,
|
||
#[serde(skip_serializing_if = "Option::is_none")]
|
||
pub duct_strength: Option<f64>,
|
||
}
|
||
|
||
#[derive(Debug, thiserror::Error)]
|
||
pub enum WriteError {
|
||
#[error("io: {0}")]
|
||
Io(#[from] std::io::Error),
|
||
#[error("encode: {0}")]
|
||
Encode(#[from] rmp_serde::encode::Error),
|
||
}
|
||
|
||
/// Build a `ScalarRow` from a single grid cell. Returns `None` for
|
||
/// cells whose surface temperature is missing or out of physically
|
||
/// plausible range — matches Elixir's `build_grid_cache_row` filter
|
||
/// (`temp_c < -80 or temp_c > 60` are dropped).
|
||
pub fn derive_row(
|
||
lat: f64,
|
||
lon: f64,
|
||
valid_time: DateTime<Utc>,
|
||
cell: &CellValues,
|
||
) -> Option<ScalarRow> {
|
||
let temp_c = cell
|
||
.get("TMP:2 m above ground")
|
||
.map(|&v| v as f64 - 273.15)?;
|
||
if !temp_c.is_finite() || !(-80.0..=60.0).contains(&temp_c) {
|
||
return None;
|
||
}
|
||
|
||
let dewpoint_c = cell.get("DPT:2 m above ground").map(|&v| v as f64 - 273.15);
|
||
let dewpoint_depression = dewpoint_c.map(|d| temp_c - d);
|
||
let surface_pressure_mb = cell.get("PRES:surface").map(|&v| v as f64 / 100.0);
|
||
let bl_height = cell.get("HPBL:surface").map(|&v| v as f64);
|
||
let pwat = cell
|
||
.get("PWAT:entire atmosphere (considered as a single layer)")
|
||
.map(|&v| v as f64);
|
||
|
||
let surface_rh = dewpoint_c.map(|d| 100.0 * sat_vap_pres(d) / sat_vap_pres(temp_c));
|
||
|
||
let levels = build_levels(cell);
|
||
let derived_min_grad = if levels.len() >= 3 {
|
||
min_refractivity_gradient(levels.clone())
|
||
} else {
|
||
None
|
||
};
|
||
let native_min_grad = cell.get("native_min_gradient").map(|&v| v as f64);
|
||
let refractivity_gradient = derived_min_grad.or(native_min_grad);
|
||
let ducting = Some(ducting_detected(refractivity_gradient));
|
||
let surface_refractivity_val = surface_refractivity(&levels);
|
||
|
||
let mut sorted: Vec<&Level> = levels.iter().collect();
|
||
sorted.sort_by(|a, b| {
|
||
b.pres_mb
|
||
.partial_cmp(&a.pres_mb)
|
||
.unwrap_or(std::cmp::Ordering::Equal)
|
||
});
|
||
|
||
let temp_850mb = level_value(&sorted, 850.0, |l| Some(l.tmpc));
|
||
let dewpoint_850mb = level_value(&sorted, 850.0, |l| l.dwpc);
|
||
let temp_700mb = level_value(&sorted, 700.0, |l| Some(l.tmpc));
|
||
let dewpoint_700mb = level_value(&sorted, 700.0, |l| l.dwpc);
|
||
|
||
let lapse_rate = compute_lapse_rate(&sorted);
|
||
let mid_lapse_rate = compute_layer_lapse_rate(&sorted, 850.0, 700.0);
|
||
let inversion_strength = compute_inversion_strength(&sorted);
|
||
let inversion_base_m = compute_inversion_base_m(&sorted);
|
||
|
||
// Ducts: the f01..f18 pipeline only stores the duct-summary scalars
|
||
// (`max_duct_thickness_m`, `duct_count`) per cell. Duct base requires
|
||
// the full Duct list, which isn't threaded through `CellValues` yet —
|
||
// matches the current Elixir behavior on Rust-produced profiles
|
||
// (their `:duct_characteristics` is also nil).
|
||
let duct_count = cell.get("duct_count").copied().unwrap_or(0.0);
|
||
let duct_strength = if duct_count > 0.0 {
|
||
cell.get("max_duct_thickness_m").map(|&v| v as f64)
|
||
} else {
|
||
None
|
||
};
|
||
|
||
Some(ScalarRow {
|
||
lat,
|
||
lon,
|
||
valid_time: valid_time.format("%Y-%m-%dT%H:%M:%SZ").to_string(),
|
||
temperature: Some(temp_c),
|
||
dewpoint_depression,
|
||
surface_rh,
|
||
surface_pressure_mb,
|
||
surface_refractivity: surface_refractivity_val,
|
||
refractivity_gradient,
|
||
bl_height,
|
||
pwat,
|
||
temp_850mb,
|
||
dewpoint_850mb,
|
||
temp_700mb,
|
||
dewpoint_700mb,
|
||
lapse_rate,
|
||
mid_lapse_rate,
|
||
inversion_strength,
|
||
inversion_base_m,
|
||
ducting,
|
||
duct_base_m: None,
|
||
duct_strength,
|
||
})
|
||
}
|
||
|
||
/// Absolute path the writer lands at for `valid_time`.
|
||
pub fn dir_for(scores_dir: &Path, valid_time: DateTime<Utc>) -> PathBuf {
|
||
let iso = valid_time.format("%Y-%m-%dT%H:%M:%SZ").to_string();
|
||
scores_dir.join(SUBDIR).join(iso)
|
||
}
|
||
|
||
/// Persist `rows` for `valid_time` as 5°×5° chunked `.mp.gz` files.
|
||
/// Mirrors the Elixir writer: the destination directory is wiped of
|
||
/// any prior chunks first so a smaller follow-up write doesn't leave
|
||
/// stale files behind. Each chunk is written `tmp + rename` for an
|
||
/// atomic NFS-friendly publish.
|
||
pub fn write_atomic(
|
||
scores_dir: &Path,
|
||
valid_time: DateTime<Utc>,
|
||
rows: &[ScalarRow],
|
||
) -> Result<PathBuf, WriteError> {
|
||
let dir = dir_for(scores_dir, valid_time);
|
||
std::fs::create_dir_all(&dir)?;
|
||
|
||
if let Ok(entries) = std::fs::read_dir(&dir) {
|
||
for entry in entries.flatten() {
|
||
let _ = std::fs::remove_file(entry.path());
|
||
}
|
||
}
|
||
|
||
let mut chunks: HashMap<(i32, i32), Vec<&ScalarRow>> = HashMap::new();
|
||
for row in rows {
|
||
let key = (chunk_band(row.lat), chunk_band(row.lon));
|
||
chunks.entry(key).or_default().push(row);
|
||
}
|
||
|
||
for ((lat_band, lon_band), chunk_rows) in chunks {
|
||
let path = dir.join(format!("{lat_band}_{lon_band}.mp.gz"));
|
||
|
||
let nanos = std::time::SystemTime::now()
|
||
.duration_since(std::time::UNIX_EPOCH)
|
||
.map(|d| d.as_nanos())
|
||
.unwrap_or(0);
|
||
let pid = std::process::id();
|
||
let mut tmp = path.clone().into_os_string();
|
||
tmp.push(format!(".tmp.{nanos}.{pid}"));
|
||
let tmp = PathBuf::from(tmp);
|
||
|
||
{
|
||
let buf = rmp_serde::to_vec_named(&chunk_rows)?;
|
||
let file = File::create(&tmp)?;
|
||
let mut gz = GzEncoder::new(BufWriter::new(file), Compression::default());
|
||
gz.write_all(&buf)?;
|
||
gz.finish()?;
|
||
}
|
||
|
||
std::fs::rename(&tmp, &path)?;
|
||
}
|
||
|
||
Ok(dir)
|
||
}
|
||
|
||
// ── Helpers ──────────────────────────────────────────────────────────
|
||
|
||
fn chunk_band(value: f64) -> i32 {
|
||
(value / CHUNK_STEP as f64).floor() as i32
|
||
}
|
||
|
||
fn build_levels(cell: &CellValues) -> Vec<Level> {
|
||
fetcher::grid_level_keys()
|
||
.iter()
|
||
.filter_map(|k| {
|
||
let t = cell.get(k.tmp.as_str()).copied()?;
|
||
let h = cell.get(k.hgt.as_str()).copied()?;
|
||
let d = cell.get(k.dpt.as_str()).copied();
|
||
Some(Level {
|
||
pres_mb: k.pres_mb,
|
||
hght_m: h as f64,
|
||
tmpc: (t as f64) - 273.15,
|
||
dwpc: d.map(|v| (v as f64) - 273.15),
|
||
})
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
/// Nearest level to `target_pres` within ±25 mb, mirroring Elixir's
|
||
/// `WeatherLayers.level_value/3`.
|
||
fn level_value(
|
||
sorted: &[&Level],
|
||
target_pres: f64,
|
||
project: impl Fn(&Level) -> Option<f64>,
|
||
) -> Option<f64> {
|
||
let nearest = sorted.iter().min_by(|a, b| {
|
||
(a.pres_mb - target_pres)
|
||
.abs()
|
||
.partial_cmp(&(b.pres_mb - target_pres).abs())
|
||
.unwrap_or(std::cmp::Ordering::Equal)
|
||
})?;
|
||
if (nearest.pres_mb - target_pres).abs() <= 25.0 {
|
||
project(nearest)
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
|
||
fn compute_lapse_rate(sorted: &[&Level]) -> Option<f64> {
|
||
if sorted.len() < 2 {
|
||
return None;
|
||
}
|
||
let surface = sorted.first()?;
|
||
let top = sorted.last()?;
|
||
let dh_km = (top.hght_m - surface.hght_m) / 1000.0;
|
||
if dh_km > 0.0 {
|
||
Some((surface.tmpc - top.tmpc) / dh_km)
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
|
||
fn compute_layer_lapse_rate(sorted: &[&Level], lower_pres: f64, upper_pres: f64) -> Option<f64> {
|
||
if lower_pres <= upper_pres {
|
||
return None;
|
||
}
|
||
let lower = sorted
|
||
.iter()
|
||
.find(|l| (l.pres_mb - lower_pres).abs() <= 25.0)?;
|
||
let upper = sorted
|
||
.iter()
|
||
.find(|l| (l.pres_mb - upper_pres).abs() <= 25.0)?;
|
||
let dh_km = (upper.hght_m - lower.hght_m) / 1000.0;
|
||
if dh_km > 0.0 {
|
||
Some((lower.tmpc - upper.tmpc) / dh_km)
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
|
||
fn compute_inversion_strength(sorted: &[&Level]) -> Option<f64> {
|
||
if sorted.is_empty() {
|
||
return Some(0.0);
|
||
}
|
||
let mut max_str = 0.0_f64;
|
||
for pair in sorted.windows(2) {
|
||
let dt = pair[1].tmpc - pair[0].tmpc;
|
||
if dt > max_str {
|
||
max_str = dt;
|
||
}
|
||
}
|
||
Some(max_str)
|
||
}
|
||
|
||
fn compute_inversion_base_m(sorted: &[&Level]) -> Option<f64> {
|
||
if sorted.is_empty() {
|
||
return None;
|
||
}
|
||
let sfc_hght = sorted.first()?.hght_m;
|
||
let mut max_str = 0.0_f64;
|
||
let mut base: Option<f64> = None;
|
||
for pair in sorted.windows(2) {
|
||
let dt = pair[1].tmpc - pair[0].tmpc;
|
||
if dt > 0.0 && dt > max_str {
|
||
max_str = dt;
|
||
base = Some(pair[0].hght_m - sfc_hght);
|
||
}
|
||
}
|
||
if max_str > 0.0 {
|
||
base
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use chrono::TimeZone;
|
||
use flate2::read::GzDecoder;
|
||
use std::io::Read;
|
||
use std::sync::Arc;
|
||
|
||
fn cell_with(items: &[(&str, f32)]) -> CellValues {
|
||
items
|
||
.iter()
|
||
.map(|&(k, v)| (Arc::<str>::from(k), v))
|
||
.collect()
|
||
}
|
||
|
||
fn surface_only_cell() -> CellValues {
|
||
cell_with(&[
|
||
("TMP:2 m above ground", 295.65), // 22.5 °C
|
||
("DPT:2 m above ground", 285.65), // 12.5 °C → depression 10 °C
|
||
("PRES:surface", 101_320.0), // 1013.2 mb
|
||
("HPBL:surface", 800.0),
|
||
(
|
||
"PWAT:entire atmosphere (considered as a single layer)",
|
||
25.0,
|
||
),
|
||
])
|
||
}
|
||
|
||
#[test]
|
||
fn drops_cells_with_unphysical_surface_temp() {
|
||
let mut cell = surface_only_cell();
|
||
cell.insert("TMP:2 m above ground".into(), 100.0); // way too cold
|
||
let vt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap();
|
||
assert!(derive_row(33.0, -97.0, vt, &cell).is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn surface_only_row_carries_basic_fields() {
|
||
let cell = surface_only_cell();
|
||
let vt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap();
|
||
let row = derive_row(33.0, -97.0, vt, &cell).expect("row");
|
||
|
||
assert_eq!(row.lat, 33.0);
|
||
assert_eq!(row.lon, -97.0);
|
||
assert_eq!(row.valid_time, "2026-04-29T12:00:00Z");
|
||
assert!((row.temperature.unwrap() - 22.5).abs() < 1e-3);
|
||
assert!((row.dewpoint_depression.unwrap() - 10.0).abs() < 1e-3);
|
||
assert!((row.surface_pressure_mb.unwrap() - 1013.2).abs() < 1e-3);
|
||
assert!(row.surface_rh.unwrap() > 30.0 && row.surface_rh.unwrap() < 70.0);
|
||
assert_eq!(row.bl_height, Some(800.0));
|
||
assert_eq!(row.pwat, Some(25.0));
|
||
|
||
// No pressure-level data → upper-air fields stay nil.
|
||
assert_eq!(row.temp_850mb, None);
|
||
assert_eq!(row.lapse_rate, None);
|
||
}
|
||
|
||
#[test]
|
||
fn upper_air_levels_derive_lapse_rate_and_850mb() {
|
||
// Need real fetcher::GRID_PRESSURE_LEVELS keys. Build directly
|
||
// from grid_level_keys() so the synthetic cell mirrors prod.
|
||
let mut cell = surface_only_cell();
|
||
for k in fetcher::grid_level_keys() {
|
||
// Plausible synthetic profile: lapse 6.5 °C/km, scale height
|
||
// ≈ 8 km. Heights derived from std atmosphere approximation.
|
||
let h_m = match k.pres_mb as i32 {
|
||
1000 => 100.0,
|
||
925 => 800.0,
|
||
850 => 1500.0,
|
||
700 => 3000.0,
|
||
500 => 5600.0,
|
||
250 => 10_400.0,
|
||
_ => 1500.0 + (1000.0 - k.pres_mb) * 8.0,
|
||
};
|
||
let t_c = 22.5 - h_m * 6.5e-3; // °C
|
||
let d_c = t_c - 5.0;
|
||
cell.insert(k.tmp.clone().into(), (t_c + 273.15) as f32);
|
||
cell.insert(k.dpt.clone().into(), (d_c + 273.15) as f32);
|
||
cell.insert(k.hgt.clone().into(), h_m as f32);
|
||
}
|
||
|
||
let vt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap();
|
||
let row = derive_row(33.0, -97.0, vt, &cell).expect("row");
|
||
|
||
let temp_850 = row.temp_850mb.expect("850 mb T");
|
||
// h_m at 850 mb is 1500 → t_c = 22.5 - 1500 * 6.5e-3 = 12.75 °C
|
||
assert!(
|
||
(temp_850 - 12.75).abs() < 0.5,
|
||
"expected ~12.75 °C at 850 mb, got {temp_850}"
|
||
);
|
||
|
||
let lapse = row.lapse_rate.expect("lapse_rate");
|
||
assert!(
|
||
(lapse - 6.5).abs() < 0.5,
|
||
"expected ~6.5 °C/km lapse rate, got {lapse}"
|
||
);
|
||
|
||
let mid = row.mid_lapse_rate.expect("mid_lapse_rate");
|
||
assert!(
|
||
(mid - 6.5).abs() < 0.5,
|
||
"expected ~6.5 mid lapse, got {mid}"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn write_atomic_round_trip_via_msgpack() {
|
||
let dir = tempfile::tempdir().unwrap();
|
||
let vt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap();
|
||
|
||
let row = ScalarRow {
|
||
lat: 33.0,
|
||
lon: -97.0,
|
||
valid_time: "2026-04-29T12:00:00Z".to_string(),
|
||
temperature: Some(22.5),
|
||
dewpoint_depression: Some(10.0),
|
||
surface_rh: Some(50.0),
|
||
ducting: Some(false),
|
||
..ScalarRow::default()
|
||
};
|
||
|
||
write_atomic(dir.path(), vt, &[row]).unwrap();
|
||
|
||
let chunk = dir
|
||
.path()
|
||
.join(SUBDIR)
|
||
.join("2026-04-29T12:00:00Z")
|
||
.join("6_-20.mp.gz");
|
||
assert!(chunk.exists());
|
||
|
||
let raw = std::fs::read(&chunk).unwrap();
|
||
let mut gz = GzDecoder::new(&raw[..]);
|
||
let mut buf = Vec::new();
|
||
gz.read_to_end(&mut buf).unwrap();
|
||
let decoded: rmpv::Value = rmp_serde::from_slice(&buf).unwrap();
|
||
let arr = decoded.as_array().expect("top-level array");
|
||
assert_eq!(arr.len(), 1);
|
||
|
||
let rmpv::Value::Map(pairs) = &arr[0] else {
|
||
panic!("row should be a map, got {:?}", arr[0]);
|
||
};
|
||
let lat = pairs
|
||
.iter()
|
||
.find(|(k, _)| k.as_str() == Some("lat"))
|
||
.unwrap();
|
||
assert_eq!(lat.1.as_f64(), Some(33.0));
|
||
|
||
let temp = pairs
|
||
.iter()
|
||
.find(|(k, _)| k.as_str() == Some("temperature"))
|
||
.unwrap();
|
||
assert_eq!(temp.1.as_f64(), Some(22.5));
|
||
}
|
||
|
||
#[test]
|
||
fn chunk_band_floors_match_elixir() {
|
||
// floor(33.0 / 5) = 6, floor(-97.0 / 5) = -20 (matches Elixir's
|
||
// `chunk_band` Float.floor/trunc).
|
||
assert_eq!(chunk_band(33.0), 6);
|
||
assert_eq!(chunk_band(-97.0), -20);
|
||
assert_eq!(chunk_band(35.0), 7);
|
||
assert_eq!(chunk_band(-95.0), -19);
|
||
}
|
||
}
|