b3ad3bdbe6
Python pskr_mqtt_listen.py: - Guard against malformed CONNACK/SUBACK/PUBLISH packets - Fix _s() truthiness: is not None instead of if v (zero SNR was lost) - EINTR-safe select loop, OOB data handling, VBInt overflow check - Proper socket cleanup with try/finally + DISCONNECT on shutdown Elixir: - Log dropped spot errors in aggregator (was silently discarded) - Wire retain_scores_window into NotifyListener chain completion - Recurse sweep_tmp_dir into band/weather_scalars subdirectories - Rescue recalibrator.run/0 to write failed status row on crash - Handle nil in fmt_snr/fmt_callsigns (no more nil dB crashes) - Replace Stream.run with Enum.reduce logging exits in poll_worker - Handle File.stat race in ms_footprints prune, grid_center nil log - Fix unused variables, stale comments, missing get_path!/1 Rust: - Graceful JoinError handling in fetcher/hrdps_fetcher (no more panics) - round_to_5min returns Option (leap-second safe) - Acquire/Release ordering on shutdown flags (was Relaxed) - Parameterized NOTIFY in db.rs, defensive scalar sweep, NaN guard - OsString::push for tmp naming, clippy fixes
1263 lines
48 KiB
Rust
1263 lines
48 KiB
Rust
//! End-to-end f01..f48 chain step: fetch surface + pressure GRIBs, decode,
|
||
//! build per-cell `Conditions`, score every band, write score-grid files
|
||
//! (`<band>/<iso>.prop`) to the shared scores directory.
|
||
//!
|
||
//! This is the Rust equivalent of `PropagationGridWorker.process_forecast_hour/4`
|
||
//! minus the f00-only enrichment (native duct, NEXRAD, commercial link
|
||
//! boost, ProfilesFile write) — Elixir retains f00.
|
||
|
||
use std::path::Path;
|
||
use std::sync::Arc;
|
||
|
||
use chrono::{DateTime, Datelike, Timelike, Utc};
|
||
|
||
use crate::band_config;
|
||
use crate::commercial::{self, LinkLookupEntry};
|
||
use crate::decoder::{self, CellValues, PointGrid};
|
||
use crate::fetcher::{self, HrrrClient, Product};
|
||
use crate::grid::{hrdps_only_points, wgrib2_grid_spec};
|
||
use crate::hrdps_fetcher::{self, HrdpsClient};
|
||
use crate::native_duct;
|
||
use crate::nexrad::{self, NexradObservation};
|
||
use crate::profiles_file::{self, CellEntry};
|
||
use crate::scorer::{self, Conditions};
|
||
use crate::scores_file::{self, ScorePoint};
|
||
use crate::sounding_params::{self, Level};
|
||
use crate::weather_scalar_file::{self, ScalarRow};
|
||
|
||
#[derive(Debug, thiserror::Error)]
|
||
pub enum PipelineError {
|
||
#[error("fetch: {0}")]
|
||
Fetch(#[from] fetcher::FetchError),
|
||
#[error("hrdps fetch: {0}")]
|
||
HrdpsFetch(#[from] hrdps_fetcher::HrdpsFetchError),
|
||
#[error("decode: {0}")]
|
||
Decode(#[from] decoder::DecodeError),
|
||
#[error("write: {0}")]
|
||
Write(#[from] scores_file::WriteError),
|
||
#[error("profile write: {0}")]
|
||
ProfileWrite(#[from] profiles_file::WriteError),
|
||
#[error("scalar write: {0}")]
|
||
ScalarWrite(#[from] weather_scalar_file::WriteError),
|
||
#[error("native duct: {0}")]
|
||
NativeDuct(#[from] native_duct::NativeDuctError),
|
||
#[error("nexrad: {0}")]
|
||
Nexrad(#[from] nexrad::NexradError),
|
||
#[error("commercial: {0}")]
|
||
Commercial(#[from] commercial::CommercialError),
|
||
#[error("forecast hour 0 is reserved for Elixir")]
|
||
F00Reserved,
|
||
#[error("surface grib was empty")]
|
||
EmptyGrib,
|
||
}
|
||
|
||
#[derive(Debug, Clone)]
|
||
pub struct ChainStepInput {
|
||
pub run_time: DateTime<Utc>,
|
||
pub forecast_hour: u8,
|
||
}
|
||
|
||
impl ChainStepInput {
|
||
pub fn valid_time(&self) -> DateTime<Utc> {
|
||
self.run_time + chrono::Duration::hours(self.forecast_hour as i64)
|
||
}
|
||
}
|
||
|
||
#[derive(Debug, Clone)]
|
||
pub struct ChainStepStats {
|
||
pub score_files_written: u32,
|
||
pub point_count: u32,
|
||
pub band_count: u32,
|
||
/// Number of cells written to the per-valid_time profile file.
|
||
/// Always `point_count` after the f01..f48 profile-file write was
|
||
/// added on 2026-04-25 — kept as a separate field so the log line
|
||
/// is symmetric with the analysis step's `profile_cells_written`
|
||
/// in `RunStepStats`, and so a future failure mode that writes
|
||
/// scores but skips the profile file shows up as a divergence.
|
||
pub profile_cells_written: u32,
|
||
}
|
||
|
||
#[tracing::instrument(
|
||
name = "pipeline.run_chain_step",
|
||
skip(client, scores_dir),
|
||
fields(
|
||
run_time = %step.run_time,
|
||
forecast_hour = step.forecast_hour,
|
||
)
|
||
)]
|
||
pub async fn run_chain_step(
|
||
client: &HrrrClient,
|
||
scores_dir: &Path,
|
||
step: &ChainStepInput,
|
||
) -> Result<ChainStepStats, PipelineError> {
|
||
if step.forecast_hour == 0 {
|
||
return Err(PipelineError::F00Reserved);
|
||
}
|
||
let valid_time = step.valid_time();
|
||
let date = step.run_time.date_naive();
|
||
let hour = step.run_time.hour() as u8;
|
||
|
||
let sfc_wanted = fetcher::surface_messages_owned();
|
||
let prs_wanted = fetcher::pressure_messages_grid();
|
||
let sfc_fut = client.fetch_product_blob(
|
||
date,
|
||
hour,
|
||
Product::Surface,
|
||
step.forecast_hour,
|
||
&sfc_wanted,
|
||
);
|
||
let prs_fut = client.fetch_product_blob(
|
||
date,
|
||
hour,
|
||
Product::Pressure,
|
||
step.forecast_hour,
|
||
&prs_wanted,
|
||
);
|
||
let (sfc_blob, prs_blob) = tokio::try_join!(sfc_fut, prs_fut)?;
|
||
if sfc_blob.is_empty() {
|
||
return Err(PipelineError::EmptyGrib);
|
||
}
|
||
|
||
let grid_spec = wgrib2_grid_spec();
|
||
let sfc_pattern = match_pattern(&sfc_wanted);
|
||
let prs_pattern = match_pattern(&prs_wanted);
|
||
|
||
// Decoder runs wgrib2 subprocess — block so we don't hog the tokio
|
||
// reactor. `spawn_blocking` lifts it onto the dedicated pool.
|
||
let grid_spec_sfc = grid_spec;
|
||
let sfc_grid = tokio::task::spawn_blocking(move || {
|
||
decoder::extract_grid(&sfc_blob, &sfc_pattern, grid_spec_sfc)
|
||
})
|
||
.await
|
||
.expect("blocking join")?;
|
||
|
||
let grid_spec_prs = grid_spec;
|
||
let prs_grid = tokio::task::spawn_blocking(move || {
|
||
decoder::extract_grid(&prs_blob, &prs_pattern, grid_spec_prs)
|
||
})
|
||
.await
|
||
.expect("blocking join")?;
|
||
|
||
let merged = merge_grids(sfc_grid, prs_grid);
|
||
let point_count = merged.len() as u32;
|
||
|
||
// Build `prepared` (for scoring) and `profile_entries` (for the
|
||
// per-valid_time profile file consumed by /skewt and the LiveView
|
||
// point-detail panels) in a single pass over `merged`. We need a
|
||
// borrowing iterator (not `into_iter`) because `cell_to_profile_entry`
|
||
// also reads from the cell — `merged` is dropped immediately after
|
||
// this loop finishes, so the ~200 MB string-keyed HashMaps don't
|
||
// outlive the scoring stage. The price is briefly retaining both
|
||
// vecs at the same time as `merged`, which is well inside the
|
||
// 512 MiB cgroup ceiling on talos5.
|
||
let mut prepared: Vec<(f64, f64, Conditions, scorer::BandInvariants)> =
|
||
Vec::with_capacity(merged.len());
|
||
let mut profile_entries: Vec<CellEntry> = Vec::with_capacity(merged.len());
|
||
let mut scalar_rows: Vec<ScalarRow> = Vec::with_capacity(merged.len());
|
||
for (key, cell) in merged.iter() {
|
||
let (lat, lon) = decoder::key_to_latlon(*key);
|
||
if let Some(conditions) = cell_to_conditions(cell, lat, lon, &valid_time) {
|
||
let invariants = scorer::precompute_band_invariants(&conditions);
|
||
prepared.push((lat, lon, conditions, invariants));
|
||
profile_entries.push(cell_to_profile_entry(lat, lon, cell));
|
||
if let Some(row) = weather_scalar_file::derive_row(lat, lon, valid_time, cell) {
|
||
scalar_rows.push(row);
|
||
}
|
||
}
|
||
}
|
||
drop(merged);
|
||
let profile_cells_written = profile_entries.len() as u32;
|
||
|
||
// Spawn the profile-file write on the blocking pool so it overlaps
|
||
// with the scoring loop instead of serialising on the hot path.
|
||
let scores_dir_owned = scores_dir.to_path_buf();
|
||
let scores_dir_for_profiles = scores_dir_owned.clone();
|
||
let profile_future = {
|
||
let profiles = profile_entries;
|
||
tokio::task::spawn_blocking(move || {
|
||
profiles_file::write_atomic(&scores_dir_for_profiles, valid_time, &profiles)
|
||
.map(|_| 0u32)
|
||
})
|
||
};
|
||
|
||
// Persist the derived scalar artifact alongside the raw profile so
|
||
// Elixir's `/weather` reads land on a pre-derived chunk file
|
||
// instead of falling back to the cold ProfilesFile decode + per-cell
|
||
// SoundingParams + WeatherLayers derivation. Identical wire format
|
||
// (gzipped MessagePack chunked by 5°×5°) to what the Elixir writer
|
||
// produces — the reader doesn't care which side wrote the bytes.
|
||
let scores_dir_for_scalars = scores_dir_owned.clone();
|
||
let scalar_future = {
|
||
let rows = scalar_rows;
|
||
tokio::task::spawn_blocking(move || {
|
||
weather_scalar_file::write_atomic(&scores_dir_for_scalars, valid_time, &rows)
|
||
.map(|_| 0u32)
|
||
})
|
||
};
|
||
|
||
// Step 2: score and write every band in parallel.
|
||
//
|
||
// Scoring: rayon's thread pool over the 23 bands. Each band is
|
||
// independent pure math across the ~92k-cell `prepared` slice;
|
||
// sharing is read-only so no locking. This cuts per-step scoring
|
||
// wall from ~20s serial to ~5-8s on a 4-core box.
|
||
//
|
||
// Writing: each band file is an independent NFS fsync. Launching
|
||
// them as parallel `spawn_blocking` tasks lets them overlap on the
|
||
// blocking pool instead of serializing one-at-a-time on the hot
|
||
// path. Memory stays bounded — each Vec<ScorePoint> is ~1.8 MB so
|
||
// 23 concurrent ≈ ~40 MB peak.
|
||
let bands = band_config::all_bands();
|
||
let prepared = Arc::new(prepared);
|
||
let scored: Vec<(u32, Vec<ScorePoint>)> = {
|
||
use rayon::prelude::*;
|
||
bands
|
||
.par_iter()
|
||
.map(|band| {
|
||
let mut scores: Vec<ScorePoint> = Vec::with_capacity(prepared.len());
|
||
for (lat, lon, conditions, invariants) in prepared.iter() {
|
||
let r = scorer::composite_score_with(conditions, band, Some(*invariants));
|
||
scores.push(ScorePoint {
|
||
lat: *lat,
|
||
lon: *lon,
|
||
score: r.score,
|
||
});
|
||
}
|
||
(band.freq_mhz, scores)
|
||
})
|
||
.collect()
|
||
};
|
||
|
||
let write_handles: Vec<_> = scored
|
||
.into_iter()
|
||
.map(|(band_mhz, scores)| {
|
||
let dir = scores_dir_owned.clone();
|
||
tokio::task::spawn_blocking(move || {
|
||
scores_file::write_atomic(&dir, band_mhz, valid_time, &scores)
|
||
})
|
||
})
|
||
.collect();
|
||
for h in write_handles {
|
||
h.await.expect("blocking join")?;
|
||
}
|
||
profile_future
|
||
.await
|
||
.expect("blocking join")
|
||
.map_err(PipelineError::ProfileWrite)?;
|
||
scalar_future
|
||
.await
|
||
.expect("blocking join")
|
||
.map_err(PipelineError::ScalarWrite)?;
|
||
let files_written = bands.len() as u32;
|
||
|
||
Ok(ChainStepStats {
|
||
score_files_written: files_written,
|
||
point_count,
|
||
band_count: bands.len() as u32,
|
||
profile_cells_written,
|
||
})
|
||
}
|
||
|
||
/// HRDPS forecast-hour chain step. Sibling of `run_chain_step` for
|
||
/// `source = 'hrdps'` rows. Differences from the HRRR path:
|
||
///
|
||
/// * Single concatenated multi-record blob from
|
||
/// `hrdps_fetcher::HrdpsClient::fetch_combined_blob` (per-variable
|
||
/// fetch + byte-concat) instead of separate sfc/prs blobs.
|
||
/// * `hrdps_grid_spec()` for the wgrib2 grid extraction (Canadian
|
||
/// bbox, rotated lat/lon handled internally by wgrib2).
|
||
/// * Cells inside HRRR's CONUS bbox are dropped from the output —
|
||
/// HRRR owns those — using a HashSet of `hrdps_only_points()` as
|
||
/// the membership test.
|
||
/// * DEPR (T-Td depression) is back-filled to DPT keys before
|
||
/// `cell_to_conditions` runs, since HRDPS publishes DEPR as a
|
||
/// primary surface variable rather than DPT.
|
||
/// * Writes to `<base>/<band>/<iso>.hrdps.prop` via
|
||
/// `scores_file::write_atomic_hrdps`. Sibling file to HRRR's
|
||
/// `.prop`, never overwrites.
|
||
/// * Skips the per-valid_time profile + scalar artifacts — those
|
||
/// would clobber the HRRR-written equivalents at the same
|
||
/// `valid_time`. Re-introducing them when HRDPS surfaces on
|
||
/// `/weather` is a separate stage.
|
||
#[tracing::instrument(
|
||
name = "pipeline.run_chain_step_hrdps",
|
||
skip(client, scores_dir),
|
||
fields(
|
||
run_time = %step.run_time,
|
||
forecast_hour = step.forecast_hour,
|
||
)
|
||
)]
|
||
pub async fn run_chain_step_hrdps(
|
||
client: &HrdpsClient,
|
||
scores_dir: &Path,
|
||
step: &ChainStepInput,
|
||
) -> Result<ChainStepStats, PipelineError> {
|
||
if step.forecast_hour == 0 {
|
||
return Err(PipelineError::F00Reserved);
|
||
}
|
||
let valid_time = step.valid_time();
|
||
let cycle = step.run_time;
|
||
|
||
let blob = client
|
||
.fetch_combined_blob(cycle, step.forecast_hour)
|
||
.await?;
|
||
if blob.is_empty() {
|
||
return Err(PipelineError::EmptyGrib);
|
||
}
|
||
|
||
// wgrib2 inventory keys for HRDPS use the same NCEP-style level
|
||
// strings as HRRR (`TMP:2 m above ground`, `DEPR:850 mb`, etc.) —
|
||
// the MSC filename convention is independent of what wgrib2 prints
|
||
// when scanning the file. Match every variable HRDPS publishes;
|
||
// post-filter happens in cell_to_conditions.
|
||
let pattern = ":(TMP|DPT|DEPR|PRES|HPBL|UGRD|VGRD|TCDC|HGT):";
|
||
|
||
// Use point-extract (`wgrib2 -lon`) instead of full-grid `-lola`
|
||
// interpolation. The full-grid path takes >10 min per chain step on
|
||
// HRDPS's rotated lat/lon source (production observation 2026-04-29);
|
||
// -lon only computes the rotation math at the points we actually
|
||
// need, dropping wall time to ~30-90 s for the ~57k Canadian cells.
|
||
let canadian_points = hrdps_only_points();
|
||
let point_count_for_log = canadian_points.len() as u32;
|
||
let pattern_owned = pattern.to_string();
|
||
let blob_for_decode = blob.clone();
|
||
let grid = tokio::task::spawn_blocking(move || {
|
||
decoder::extract_points(&blob_for_decode, &pattern_owned, &canadian_points)
|
||
})
|
||
.await
|
||
.expect("blocking join")?;
|
||
drop(blob);
|
||
|
||
let point_count = grid.len() as u32;
|
||
tracing::info!(
|
||
requested = point_count_for_log,
|
||
decoded = point_count,
|
||
"hrdps points extracted"
|
||
);
|
||
|
||
let mut prepared: Vec<(f64, f64, Conditions, scorer::BandInvariants)> =
|
||
Vec::with_capacity(grid.len());
|
||
let mut scalar_rows: Vec<ScalarRow> = Vec::with_capacity(grid.len());
|
||
for (key, mut cell) in grid.into_iter() {
|
||
let (lat, lon) = decoder::key_to_latlon(key);
|
||
// HRDPS publishes DEPR (T-Td depression in K) as a primary
|
||
// surface + pressure-level variable rather than DPT. Back-fill
|
||
// DPT entries from `TMP - DEPR` so cell_to_conditions sees the
|
||
// expected key set.
|
||
backfill_dpt_from_depr(&mut cell);
|
||
if let Some(conditions) = cell_to_conditions(&cell, lat, lon, &valid_time) {
|
||
let invariants = scorer::precompute_band_invariants(&conditions);
|
||
prepared.push((lat, lon, conditions, invariants));
|
||
// Derive the per-cell scalar row alongside scoring. The /weather
|
||
// map reads these via Weather.weather_grid_at to render
|
||
// temperature/dewpoint-depression/refractivity layers — surfacing
|
||
// HRDPS here is what gives Canadian cells weather data on the
|
||
// map north of the HRRR coverage line.
|
||
if let Some(row) = weather_scalar_file::derive_row(lat, lon, valid_time, &cell) {
|
||
scalar_rows.push(row);
|
||
}
|
||
}
|
||
}
|
||
|
||
let bands = band_config::all_bands();
|
||
let prepared = std::sync::Arc::new(prepared);
|
||
let scored: Vec<(u32, Vec<ScorePoint>)> = {
|
||
use rayon::prelude::*;
|
||
bands
|
||
.par_iter()
|
||
.map(|band| {
|
||
let mut scores: Vec<ScorePoint> = Vec::with_capacity(prepared.len());
|
||
for (lat, lon, conditions, invariants) in prepared.iter() {
|
||
let r = scorer::composite_score_with(conditions, band, Some(*invariants));
|
||
scores.push(ScorePoint {
|
||
lat: *lat,
|
||
lon: *lon,
|
||
score: r.score,
|
||
});
|
||
}
|
||
(band.freq_mhz, scores)
|
||
})
|
||
.collect()
|
||
};
|
||
|
||
let scores_dir_owned = scores_dir.to_path_buf();
|
||
|
||
// Write the HRDPS scalar artifact alongside the scores. Sibling dir
|
||
// (`<vt>.hrdps/`) so HRRR's `<vt>/` write isn't clobbered. Without
|
||
// this the /weather map would only show CONUS data.
|
||
let scalar_dir = scores_dir_owned.clone();
|
||
let scalar_future = tokio::task::spawn_blocking(move || {
|
||
weather_scalar_file::write_atomic_hrdps(&scalar_dir, valid_time, &scalar_rows).map(|_| ())
|
||
});
|
||
|
||
let write_handles: Vec<_> = scored
|
||
.into_iter()
|
||
.map(|(band_mhz, scores)| {
|
||
let dir = scores_dir_owned.clone();
|
||
tokio::task::spawn_blocking(move || {
|
||
scores_file::write_atomic_hrdps(&dir, band_mhz, valid_time, &scores)
|
||
})
|
||
})
|
||
.collect();
|
||
for h in write_handles {
|
||
h.await.expect("blocking join")?;
|
||
}
|
||
scalar_future
|
||
.await
|
||
.expect("blocking join")
|
||
.map_err(PipelineError::ScalarWrite)?;
|
||
let files_written = bands.len() as u32;
|
||
|
||
Ok(ChainStepStats {
|
||
score_files_written: files_written,
|
||
point_count,
|
||
band_count: bands.len() as u32,
|
||
profile_cells_written: 0,
|
||
})
|
||
}
|
||
|
||
fn backfill_dpt_from_depr(cell: &mut CellValues) {
|
||
// Snapshot of (TMP key, DEPR key, DPT key) tuples to insert. Done in
|
||
// two passes so the iteration in pass 1 isn't invalidated by the
|
||
// mutation in pass 2.
|
||
let pending: Vec<(std::sync::Arc<str>, f32)> = cell
|
||
.keys()
|
||
.filter_map(|key| {
|
||
let depr_key = key.strip_prefix("DEPR:")?;
|
||
let tmp_key: std::sync::Arc<str> = format!("TMP:{depr_key}").into();
|
||
let dpt_key: std::sync::Arc<str> = format!("DPT:{depr_key}").into();
|
||
if cell.contains_key(&dpt_key) {
|
||
return None;
|
||
}
|
||
let tmp = cell.get(&tmp_key).copied()?;
|
||
let depr = cell.get(key).copied()?;
|
||
Some((dpt_key, tmp - depr))
|
||
})
|
||
.collect();
|
||
|
||
for (dpt_key, dpt_val) in pending {
|
||
cell.insert(dpt_key, dpt_val);
|
||
}
|
||
}
|
||
|
||
/// wgrib2 `-match` regex: `":(A|B|C):"`. Var names only — levels are
|
||
/// post-filtered when we read the binary back.
|
||
fn match_pattern(wanted: &[(String, String)]) -> String {
|
||
let mut vars: Vec<&str> = wanted.iter().map(|(v, _)| v.as_str()).collect();
|
||
vars.sort();
|
||
vars.dedup();
|
||
format!(":({}):", vars.join("|"))
|
||
}
|
||
|
||
fn merge_grids(mut sfc: PointGrid, prs: PointGrid) -> PointGrid {
|
||
for (k, v) in prs {
|
||
sfc.entry(k).or_default().extend(v);
|
||
}
|
||
sfc
|
||
}
|
||
|
||
fn cell_to_conditions(
|
||
cell: &CellValues,
|
||
lat: f64,
|
||
lon: f64,
|
||
valid_time: &DateTime<Utc>,
|
||
) -> Option<Conditions> {
|
||
let tmp_k = cell.get("TMP:2 m above ground").copied()?;
|
||
let dpt_k = cell.get("DPT:2 m above ground").copied()?;
|
||
let temp_c = (tmp_k - 273.15) as f64;
|
||
let dewpoint_c = (dpt_k - 273.15) as f64;
|
||
let temp_f = scorer::c_to_f(temp_c);
|
||
let dewpoint_f = scorer::c_to_f(dewpoint_c);
|
||
let abs_humidity = scorer::absolute_humidity(temp_c, dewpoint_c);
|
||
|
||
let wind_speed_kts = match (
|
||
cell.get("UGRD:10 m above ground").copied(),
|
||
cell.get("VGRD:10 m above ground").copied(),
|
||
) {
|
||
(Some(u), Some(v)) => Some(scorer::wind_speed_kts(u as f64, v as f64)),
|
||
_ => None,
|
||
};
|
||
|
||
let sky_cover_pct = cell
|
||
.get("TCDC:entire atmosphere")
|
||
.copied()
|
||
.map(|v| v as f64);
|
||
let pwat_mm = cell
|
||
.get("PWAT:entire atmosphere (considered as a single layer)")
|
||
.copied()
|
||
.map(|v| v as f64);
|
||
let pressure_mb = cell
|
||
.get("PRES:surface")
|
||
.copied()
|
||
.map(|pa| pa as f64 / 100.0);
|
||
// Mirror Elixir's merged_rain_rate: pick the heavier of HRRR's
|
||
// accumulation-derived rate and NEXRAD's reflectivity-derived rate
|
||
// so a fast convective cell that hasn't yet shown up in the hourly
|
||
// APCP still triggers the rain penalty.
|
||
let hrrr_rate = cell
|
||
.get("APCP:surface")
|
||
.copied()
|
||
.map(|mm| mm as f64)
|
||
.unwrap_or(0.0);
|
||
let nexrad_rate = scorer::dbz_to_rain_rate_mmhr(
|
||
cell.get("nexrad_max_reflectivity_dbz")
|
||
.copied()
|
||
.map(|v| v as f64),
|
||
);
|
||
let merged_rate = hrrr_rate.max(nexrad_rate);
|
||
let rain_rate_mmhr = if merged_rate > 0.0 {
|
||
Some(merged_rate)
|
||
} else {
|
||
None
|
||
};
|
||
let bl_depth_m = cell.get("HPBL:surface").copied().map(|v| v as f64);
|
||
let best_duct_band_ghz = cell.get("best_duct_freq_ghz").copied().map(|v| v as f64);
|
||
|
||
let levels: Vec<Level> = fetcher::grid_level_keys()
|
||
.iter()
|
||
.filter_map(|k| {
|
||
let t = cell.get(k.tmp.as_str()).copied()?;
|
||
let d = cell.get(k.dpt.as_str()).copied();
|
||
let h = cell.get(k.hgt.as_str()).copied()?;
|
||
Some(Level {
|
||
pres_mb: k.pres_mb,
|
||
hght_m: h as f64,
|
||
tmpc: (t - 273.15) as f64,
|
||
dwpc: d.map(|v| (v - 273.15) as f64),
|
||
})
|
||
})
|
||
.collect();
|
||
let min_refractivity_gradient = sounding_params::min_refractivity_gradient(levels);
|
||
|
||
Some(Conditions {
|
||
abs_humidity,
|
||
temp_f,
|
||
dewpoint_f,
|
||
wind_speed_kts,
|
||
sky_cover_pct,
|
||
utc_hour: valid_time.hour() as u8,
|
||
utc_minute: valid_time.minute() as u8,
|
||
month: valid_time.month() as u8,
|
||
longitude: lon,
|
||
latitude: Some(lat),
|
||
pressure_mb,
|
||
prev_pressure_mb: None,
|
||
rain_rate_mmhr,
|
||
min_refractivity_gradient,
|
||
bl_depth_m,
|
||
pwat_mm,
|
||
best_duct_band_ghz,
|
||
bulk_richardson: None,
|
||
})
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
#[test]
|
||
fn match_pattern_dedups_vars() {
|
||
let wanted = vec![
|
||
("TMP".into(), "2 m above ground".into()),
|
||
("DPT".into(), "2 m above ground".into()),
|
||
("TMP".into(), "1000 mb".into()),
|
||
];
|
||
assert_eq!(match_pattern(&wanted), ":(DPT|TMP):");
|
||
}
|
||
|
||
#[test]
|
||
fn cell_to_conditions_maps_fields() {
|
||
use chrono::TimeZone;
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 298.15); // 25 °C
|
||
cell.insert("DPT:2 m above ground".into(), 293.15); // 20 °C
|
||
cell.insert("UGRD:10 m above ground".into(), 3.0);
|
||
cell.insert("VGRD:10 m above ground".into(), 4.0);
|
||
cell.insert("TCDC:entire atmosphere".into(), 30.0);
|
||
cell.insert("PRES:surface".into(), 101_000.0);
|
||
cell.insert("HPBL:surface".into(), 400.0);
|
||
cell.insert(
|
||
"PWAT:entire atmosphere (considered as a single layer)".into(),
|
||
25.0,
|
||
);
|
||
|
||
let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
|
||
let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
|
||
assert!((c.temp_f - 77.0).abs() < 0.1);
|
||
assert!((c.dewpoint_f - 68.0).abs() < 0.1);
|
||
assert!((c.wind_speed_kts.unwrap() - 9.72).abs() < 0.1);
|
||
assert!((c.pressure_mb.unwrap() - 1010.0).abs() < 0.01);
|
||
assert_eq!(c.sky_cover_pct, Some(30.0));
|
||
assert_eq!(c.bl_depth_m, Some(400.0));
|
||
assert_eq!(c.pwat_mm, Some(25.0));
|
||
assert_eq!(c.month, 6);
|
||
assert_eq!(c.longitude, -97.0);
|
||
}
|
||
|
||
#[test]
|
||
fn cell_missing_surface_returns_none() {
|
||
use chrono::TimeZone;
|
||
let cell = CellValues::new();
|
||
let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
|
||
assert!(cell_to_conditions(&cell, 32.0, -97.0, &vt).is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn nexrad_reflectivity_lifts_rain_rate_above_hrrr_precip() {
|
||
// Mirrors the Elixir scorer's rain merge: when HRRR's hourly
|
||
// accumulation hasn't caught a fast convective cell but NEXRAD
|
||
// sees 45 dBZ overhead, take whichever rate is higher. Without
|
||
// this the Rust pipeline silently underestimates rain attenuation.
|
||
use chrono::TimeZone;
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 295.0);
|
||
cell.insert("DPT:2 m above ground".into(), 285.0);
|
||
cell.insert("APCP:surface".into(), 0.0);
|
||
cell.insert("nexrad_max_reflectivity_dbz".into(), 45.0);
|
||
|
||
let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
|
||
let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
|
||
let rate = c.rain_rate_mmhr.expect("nexrad-derived rain rate present");
|
||
assert!(rate > 5.0, "expected >5 mm/hr from 45 dBZ, got {rate}");
|
||
}
|
||
|
||
#[test]
|
||
fn duct_freq_threads_into_best_duct_band_ghz() {
|
||
// After native_duct::merge_duct_grid runs, every cell with a
|
||
// detected trapping layer carries `best_duct_freq_ghz`. The
|
||
// scorer reads `best_duct_band_ghz` from Conditions for the 15%
|
||
// Native Duct Boost. cell_to_conditions must wire one to the
|
||
// other or the boost never fires from the Rust pipeline.
|
||
use chrono::TimeZone;
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 295.0);
|
||
cell.insert("DPT:2 m above ground".into(), 285.0);
|
||
cell.insert("best_duct_freq_ghz".into(), 24.0);
|
||
|
||
let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
|
||
let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
|
||
assert_eq!(c.best_duct_band_ghz, Some(24.0));
|
||
}
|
||
|
||
#[test]
|
||
fn rejects_forecast_hour_zero() {
|
||
use chrono::TimeZone;
|
||
let client = HrrrClient::default_base().unwrap();
|
||
let dir = tempfile::tempdir().unwrap();
|
||
let step = ChainStepInput {
|
||
run_time: Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(),
|
||
forecast_hour: 0,
|
||
};
|
||
let rt = tokio::runtime::Runtime::new().unwrap();
|
||
let err = rt
|
||
.block_on(run_chain_step(&client, dir.path(), &step))
|
||
.unwrap_err();
|
||
assert!(matches!(err, PipelineError::F00Reserved));
|
||
}
|
||
|
||
#[test]
|
||
fn hrdps_rejects_forecast_hour_zero() {
|
||
use chrono::TimeZone;
|
||
let client = HrdpsClient::default_base().unwrap();
|
||
let dir = tempfile::tempdir().unwrap();
|
||
let step = ChainStepInput {
|
||
run_time: Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(),
|
||
forecast_hour: 0,
|
||
};
|
||
let rt = tokio::runtime::Runtime::new().unwrap();
|
||
let err = rt
|
||
.block_on(run_chain_step_hrdps(&client, dir.path(), &step))
|
||
.unwrap_err();
|
||
assert!(matches!(err, PipelineError::F00Reserved));
|
||
}
|
||
|
||
#[test]
|
||
fn backfill_dpt_from_depr_derives_missing_dewpoint() {
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 290.15); // 17°C
|
||
cell.insert("DEPR:2 m above ground".into(), 5.0);
|
||
|
||
backfill_dpt_from_depr(&mut cell);
|
||
|
||
let dpt = cell
|
||
.get("DPT:2 m above ground")
|
||
.copied()
|
||
.expect("DPT back-filled from DEPR");
|
||
assert!((dpt - 285.15).abs() < 1e-3, "got {dpt}");
|
||
}
|
||
|
||
#[test]
|
||
fn backfill_dpt_does_not_clobber_existing_dpt() {
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 290.15);
|
||
cell.insert("DEPR:2 m above ground".into(), 5.0);
|
||
cell.insert("DPT:2 m above ground".into(), 280.0); // pre-existing
|
||
|
||
backfill_dpt_from_depr(&mut cell);
|
||
|
||
assert_eq!(cell.get("DPT:2 m above ground").copied(), Some(280.0));
|
||
}
|
||
|
||
#[test]
|
||
fn backfill_dpt_handles_pressure_levels() {
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:850 mb".into(), 275.0);
|
||
cell.insert("DEPR:850 mb".into(), 4.0);
|
||
cell.insert("TMP:700 mb".into(), 265.0);
|
||
cell.insert("DEPR:700 mb".into(), 3.0);
|
||
|
||
backfill_dpt_from_depr(&mut cell);
|
||
|
||
assert_eq!(cell.get("DPT:850 mb").copied(), Some(271.0));
|
||
assert_eq!(cell.get("DPT:700 mb").copied(), Some(262.0));
|
||
}
|
||
|
||
/// Integration-ish: hand-build a 2-cell surface + pressure grid,
|
||
/// merge → cell_to_conditions → scorer → scores_file::write_atomic
|
||
/// → scores_file::read, asserting the full chain agrees with
|
||
/// itself. Covers the post-fetch hot path (the one that scales
|
||
/// with grid size) without standing up wgrib2.
|
||
#[test]
|
||
fn merge_score_write_read_roundtrip() {
|
||
use chrono::TimeZone;
|
||
|
||
// decoder encodes lat/lon at millidegree precision in the key.
|
||
fn latlon_to_key(lat: f64, lon: f64) -> (i32, i32) {
|
||
((lat * 1000.0).round() as i32, (lon * 1000.0).round() as i32)
|
||
}
|
||
|
||
let mut sfc = PointGrid::new();
|
||
let mut prs = PointGrid::new();
|
||
|
||
let cells = [(32.0_f64, -97.0_f64), (32.03_f64, -97.03_f64)];
|
||
for &(lat, lon) in &cells {
|
||
let mut s = CellValues::new();
|
||
s.insert("TMP:2 m above ground".into(), 295.0);
|
||
s.insert("DPT:2 m above ground".into(), 285.0);
|
||
s.insert("UGRD:10 m above ground".into(), 2.0);
|
||
s.insert("VGRD:10 m above ground".into(), 1.0);
|
||
s.insert("TCDC:entire atmosphere".into(), 20.0);
|
||
s.insert("PRES:surface".into(), 101_000.0);
|
||
s.insert("HPBL:surface".into(), 500.0);
|
||
s.insert(
|
||
"PWAT:entire atmosphere (considered as a single layer)".into(),
|
||
20.0,
|
||
);
|
||
sfc.insert(latlon_to_key(lat, lon), s);
|
||
|
||
let mut p = CellValues::new();
|
||
p.insert("TMP:1000 mb".into(), 295.0);
|
||
p.insert("DPT:1000 mb".into(), 285.0);
|
||
p.insert("HGT:1000 mb".into(), 100.0);
|
||
p.insert("TMP:975 mb".into(), 293.0);
|
||
p.insert("DPT:975 mb".into(), 280.0);
|
||
p.insert("HGT:975 mb".into(), 300.0);
|
||
prs.insert(latlon_to_key(lat, lon), p);
|
||
}
|
||
|
||
let merged = merge_grids(sfc, prs);
|
||
assert_eq!(merged.len(), 2);
|
||
|
||
let valid_time = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
|
||
|
||
let prepared: Vec<_> = merged
|
||
.into_iter()
|
||
.filter_map(|(key, cell)| {
|
||
let (lat, lon) = crate::decoder::key_to_latlon(key);
|
||
let c = cell_to_conditions(&cell, lat, lon, &valid_time)?;
|
||
let inv = scorer::precompute_band_invariants(&c);
|
||
Some((lat, lon, c, inv))
|
||
})
|
||
.collect();
|
||
assert_eq!(prepared.len(), 2);
|
||
|
||
let band = &band_config::all_bands()[0];
|
||
let scores: Vec<ScorePoint> = prepared
|
||
.iter()
|
||
.map(|(lat, lon, c, inv)| {
|
||
let r = scorer::composite_score_with(c, band, Some(*inv));
|
||
ScorePoint {
|
||
lat: *lat,
|
||
lon: *lon,
|
||
score: r.score,
|
||
}
|
||
})
|
||
.collect();
|
||
assert_eq!(scores.len(), 2);
|
||
for s in &scores {
|
||
assert!(s.score <= 100);
|
||
}
|
||
|
||
// Write to an atomic .prop file and read it back.
|
||
let dir = tempfile::tempdir().unwrap();
|
||
scores_file::write_atomic(dir.path(), band.freq_mhz, valid_time, &scores).unwrap();
|
||
|
||
let path = scores_file::path_for(dir.path(), band.freq_mhz, valid_time);
|
||
let bytes = std::fs::read(&path).expect("read back score file");
|
||
let decoded = scores_file::decode(&bytes).expect("decode");
|
||
// Decoded grid covers the full CONUS bounding box as a dense
|
||
// byte array; our 2 score points land at their respective
|
||
// cells and the rest are sentinel no-data. Spot-check: header
|
||
// fields survived the round trip and the body is the expected
|
||
// size for a row_major n_rows × n_cols grid.
|
||
assert_eq!(decoded.band_mhz, band.freq_mhz);
|
||
assert_eq!(
|
||
decoded.body.len(),
|
||
decoded.n_rows as usize * decoded.n_cols as usize
|
||
);
|
||
}
|
||
|
||
/// cell_to_profile_entry includes pre-computed surface_refractivity
|
||
/// and min_refractivity_gradient so Elixir consumers (Skew-T,
|
||
/// contact detail) have a guaranteed fallback when SoundingParams.derive
|
||
/// hits a level without dewpoint.
|
||
#[test]
|
||
fn profile_entry_includes_refractivity_scalars() {
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 298.15);
|
||
cell.insert("DPT:2 m above ground".into(), 293.15);
|
||
cell.insert("PRES:surface".into(), 101_300.0);
|
||
cell.insert("HPBL:surface".into(), 1500.0);
|
||
cell.insert(
|
||
"PWAT:entire atmosphere (considered as a single layer)".into(),
|
||
30.0,
|
||
);
|
||
cell.insert("TMP:1000 mb".into(), 295.0);
|
||
cell.insert("DPT:1000 mb".into(), 290.0);
|
||
cell.insert("HGT:1000 mb".into(), 100.0);
|
||
cell.insert("TMP:925 mb".into(), 292.0);
|
||
cell.insert("DPT:925 mb".into(), 285.0);
|
||
cell.insert("HGT:925 mb".into(), 800.0);
|
||
cell.insert("TMP:850 mb".into(), 289.0);
|
||
cell.insert("DPT:850 mb".into(), 280.0);
|
||
cell.insert("HGT:850 mb".into(), 1500.0);
|
||
|
||
let entry = cell_to_profile_entry(32.9, -97.0, &cell);
|
||
let map = entry.profile.as_map().expect("profile is a map");
|
||
let has_sr = map.iter().any(|(k, _)| k.as_str() == Some("surface_refractivity"));
|
||
let has_mg = map.iter().any(|(k, _)| k.as_str() == Some("min_refractivity_gradient"));
|
||
assert!(has_sr, "surface_refractivity key missing");
|
||
assert!(has_mg, "min_refractivity_gradient key missing");
|
||
}
|
||
|
||
/// When every pressure level is missing dewpoint, surface_refractivity
|
||
/// and min_refractivity_gradient are omitted.
|
||
#[test]
|
||
fn profile_entry_omits_refractivity_when_no_dewpoint() {
|
||
let mut cell = CellValues::new();
|
||
cell.insert("TMP:2 m above ground".into(), 298.15);
|
||
cell.insert("PRES:surface".into(), 101_300.0);
|
||
cell.insert("TMP:1000 mb".into(), 295.0);
|
||
cell.insert("HGT:1000 mb".into(), 100.0);
|
||
cell.insert("TMP:925 mb".into(), 292.0);
|
||
cell.insert("HGT:925 mb".into(), 800.0);
|
||
|
||
let entry = cell_to_profile_entry(32.9, -97.0, &cell);
|
||
let map = entry.profile.as_map().expect("profile is a map");
|
||
let has_sr = map.iter().any(|(k, _)| k.as_str() == Some("surface_refractivity"));
|
||
let has_mg = map.iter().any(|(k, _)| k.as_str() == Some("min_refractivity_gradient"));
|
||
assert!(!has_sr, "surface_refractivity should be absent without dewpoint");
|
||
assert!(!has_mg, "min_refractivity_gradient should be absent without dewpoint");
|
||
}
|
||
}
|
||
|
||
// =====================================================================
|
||
// Analysis pipeline (f00): fetch surface+pressure + native duct, merge
|
||
// with NEXRAD + commercial-link degradation, score all bands, write the
|
||
// ProfilesFile (MessagePack) and the band score files.
|
||
// =====================================================================
|
||
|
||
#[derive(Debug, Clone)]
|
||
pub struct AnalysisStepInput {
|
||
pub run_time: DateTime<Utc>,
|
||
}
|
||
|
||
impl AnalysisStepInput {
|
||
pub fn valid_time(&self) -> DateTime<Utc> {
|
||
// f00 analysis valid_time == run_time.
|
||
self.run_time
|
||
}
|
||
}
|
||
|
||
#[derive(Debug, Clone)]
|
||
pub struct AnalysisStepStats {
|
||
pub score_files_written: u32,
|
||
pub point_count: u32,
|
||
pub band_count: u32,
|
||
pub profile_cells_written: u32,
|
||
pub duct_cells: u32,
|
||
pub nexrad_cells_with_echo: u32,
|
||
pub commercial_cells_boosted: u32,
|
||
}
|
||
|
||
#[tracing::instrument(
|
||
name = "pipeline.run_analysis_step",
|
||
skip(client, pool, scores_dir),
|
||
fields(run_time = %step.run_time)
|
||
)]
|
||
pub async fn run_analysis_step(
|
||
client: &HrrrClient,
|
||
pool: &sqlx::PgPool,
|
||
scores_dir: &Path,
|
||
step: &AnalysisStepInput,
|
||
) -> Result<AnalysisStepStats, PipelineError> {
|
||
let valid_time = step.valid_time();
|
||
let date = step.run_time.date_naive();
|
||
let hour = step.run_time.hour() as u8;
|
||
let grid_spec = wgrib2_grid_spec();
|
||
|
||
// Parallel fetch + decode of surface + pressure + native duct. All
|
||
// three are independent HRRR products; overlap the HTTP legs and
|
||
// the subsequent wgrib2 decodes so the long pole is a single
|
||
// decode, not three serial decodes.
|
||
let sfc_wanted = fetcher::surface_messages_owned();
|
||
let prs_wanted = fetcher::pressure_messages_grid();
|
||
let sfc_pattern = match_pattern(&sfc_wanted);
|
||
let prs_pattern = match_pattern(&prs_wanted);
|
||
|
||
// All three futures return different error types; lift into
|
||
// PipelineError inside each branch so try_join! stays happy.
|
||
let sfc_fut = async {
|
||
client
|
||
.fetch_product_blob(date, hour, Product::Surface, 0, &sfc_wanted)
|
||
.await
|
||
.map_err(PipelineError::Fetch)
|
||
};
|
||
let prs_fut = async {
|
||
client
|
||
.fetch_product_blob(date, hour, Product::Pressure, 0, &prs_wanted)
|
||
.await
|
||
.map_err(PipelineError::Fetch)
|
||
};
|
||
let duct_fut = async {
|
||
native_duct::fetch_native_duct_grid(client, date, hour, 0, grid_spec)
|
||
.await
|
||
.map_err(PipelineError::NativeDuct)
|
||
};
|
||
|
||
let (sfc_blob, prs_blob, duct_map) = tokio::try_join!(sfc_fut, prs_fut, duct_fut)?;
|
||
|
||
let duct_count = duct_map.len() as u32;
|
||
|
||
let sfc_grid = tokio::task::spawn_blocking(move || {
|
||
decoder::extract_grid(&sfc_blob, &sfc_pattern, grid_spec)
|
||
})
|
||
.await
|
||
.expect("blocking join")?;
|
||
|
||
let prs_grid = tokio::task::spawn_blocking(move || {
|
||
decoder::extract_grid(&prs_blob, &prs_pattern, grid_spec)
|
||
})
|
||
.await
|
||
.expect("blocking join")?;
|
||
|
||
// Merge pressure into surface, then fold duct metrics into every
|
||
// cell that has both. The duct values live alongside the raw
|
||
// CellValues strings so build_grid_cache_rows can pick them up
|
||
// under their atom keys after the Elixir reader normalises.
|
||
let mut merged = merge_grids(sfc_grid, prs_grid);
|
||
merged = native_duct::merge_duct_grid(merged, &duct_map);
|
||
|
||
// NEXRAD composite reflectivity overlay — only valid for f00
|
||
// because forecast hours can't see the future radar image.
|
||
let nexrad_http = reqwest::Client::builder()
|
||
.user_agent("prop-grid-rs/0.1")
|
||
.build()
|
||
.unwrap_or_else(|e| {
|
||
tracing::warn!(
|
||
"failed to build Nexrad HTTP client (missing root CA certs?): {} — \
|
||
f00 overlay will be skipped for this cycle",
|
||
e
|
||
);
|
||
// Return a client that will fail every request rather than
|
||
// panicking the worker. The f00 overlay is non-critical.
|
||
reqwest::Client::new()
|
||
});
|
||
let nexrad_points: Vec<(f64, f64)> =
|
||
merged.keys().map(|&k| decoder::key_to_latlon(k)).collect();
|
||
let nexrad_obs: Vec<NexradObservation> = match nexrad::fetch_frame(
|
||
&nexrad_http,
|
||
valid_time,
|
||
&nexrad_points,
|
||
)
|
||
.await
|
||
{
|
||
Ok(obs) => obs,
|
||
Err(e) => {
|
||
tracing::warn!(error = %e, "NEXRAD fetch failed — continuing without radar overlay");
|
||
Vec::new()
|
||
}
|
||
};
|
||
let nexrad_cells_with_echo = apply_nexrad(&mut merged, &nexrad_obs);
|
||
|
||
// Commercial-link degradation — precompute once, then fold into
|
||
// every in-range cell. Only ~7 links cluster around DFW so most
|
||
// cells stay untouched.
|
||
let commercial_lookup: Vec<LinkLookupEntry> = commercial::build_link_lookup(pool, valid_time)
|
||
.await
|
||
.map_err(PipelineError::Commercial)?;
|
||
let commercial_cells_boosted = apply_commercial(&mut merged, &commercial_lookup);
|
||
|
||
let point_count = merged.len() as u32;
|
||
|
||
// Build the ProfilesFile cell list (the raw enriched grid) and
|
||
// the Conditions vector (for scoring) in one pass. Profile keeps
|
||
// atmospheric + duct + NEXRAD + commercial fields under
|
||
// string-keyed msgpack; Rust doesn't need to mutate those further.
|
||
let mut profile_entries: Vec<CellEntry> = Vec::with_capacity(merged.len());
|
||
let mut prepared: Vec<(f64, f64, Conditions, scorer::BandInvariants)> =
|
||
Vec::with_capacity(merged.len());
|
||
let mut scalar_rows: Vec<ScalarRow> = Vec::with_capacity(merged.len());
|
||
for (key, cell) in merged.iter() {
|
||
let (lat, lon) = decoder::key_to_latlon(*key);
|
||
if let Some(conditions) = cell_to_conditions(cell, lat, lon, &valid_time) {
|
||
let invariants = scorer::precompute_band_invariants(&conditions);
|
||
prepared.push((lat, lon, conditions, invariants));
|
||
profile_entries.push(cell_to_profile_entry(lat, lon, cell));
|
||
if let Some(row) = weather_scalar_file::derive_row(lat, lon, valid_time, cell) {
|
||
scalar_rows.push(row);
|
||
}
|
||
}
|
||
}
|
||
|
||
// Write the profile file on the blocking pool — msgpack encode +
|
||
// gzip + atomic rename. One file per run, ~2 MB compressed, so the
|
||
// write runs well inside the scoring loop's wall time.
|
||
let scores_dir_owned = scores_dir.to_path_buf();
|
||
let scores_dir_for_profiles = scores_dir_owned.clone();
|
||
let profile_future = {
|
||
let profiles = profile_entries;
|
||
tokio::task::spawn_blocking(move || {
|
||
profiles_file::write_atomic(&scores_dir_for_profiles, valid_time, &profiles)
|
||
.map(|_| 0u32)
|
||
})
|
||
};
|
||
|
||
// Same wire format as Elixir's `ScalarFile`. See run_chain_step
|
||
// (forecast hours) for the full rationale — both code paths go
|
||
// through the same writer so f00 and f01..f48 land identical
|
||
// chunked artifacts on NFS.
|
||
let scores_dir_for_scalars = scores_dir_owned.clone();
|
||
let scalar_future = {
|
||
let rows = scalar_rows;
|
||
tokio::task::spawn_blocking(move || {
|
||
weather_scalar_file::write_atomic(&scores_dir_for_scalars, valid_time, &rows)
|
||
.map(|_| 0u32)
|
||
})
|
||
};
|
||
|
||
// Bands and score writes run in parallel, same shape as
|
||
// run_chain_step. Rayon across bands, try_join_all across writes.
|
||
let bands = band_config::all_bands();
|
||
let prepared_arc = Arc::new(prepared);
|
||
let scored: Vec<(u32, Vec<ScorePoint>)> = {
|
||
use rayon::prelude::*;
|
||
bands
|
||
.par_iter()
|
||
.map(|band| {
|
||
let mut scores: Vec<ScorePoint> = Vec::with_capacity(prepared_arc.len());
|
||
for (lat, lon, conditions, invariants) in prepared_arc.iter() {
|
||
let r = scorer::composite_score_with(conditions, band, Some(*invariants));
|
||
scores.push(ScorePoint {
|
||
lat: *lat,
|
||
lon: *lon,
|
||
score: r.score,
|
||
});
|
||
}
|
||
(band.freq_mhz, scores)
|
||
})
|
||
.collect()
|
||
};
|
||
|
||
let write_handles: Vec<_> = scored
|
||
.into_iter()
|
||
.map(|(band_mhz, scores)| {
|
||
let dir = scores_dir_owned.clone();
|
||
tokio::task::spawn_blocking(move || {
|
||
scores_file::write_atomic(&dir, band_mhz, valid_time, &scores)
|
||
})
|
||
})
|
||
.collect();
|
||
for h in write_handles {
|
||
h.await.expect("blocking join")?;
|
||
}
|
||
profile_future
|
||
.await
|
||
.expect("blocking join")
|
||
.map_err(PipelineError::ProfileWrite)?;
|
||
scalar_future
|
||
.await
|
||
.expect("blocking join")
|
||
.map_err(PipelineError::ScalarWrite)?;
|
||
|
||
let files_written = bands.len() as u32;
|
||
Ok(AnalysisStepStats {
|
||
score_files_written: files_written,
|
||
point_count,
|
||
band_count: bands.len() as u32,
|
||
profile_cells_written: prepared_arc.len() as u32,
|
||
duct_cells: duct_count,
|
||
nexrad_cells_with_echo,
|
||
commercial_cells_boosted,
|
||
})
|
||
}
|
||
|
||
fn apply_nexrad(grid: &mut PointGrid, obs: &[NexradObservation]) -> u32 {
|
||
let mut with_echo = 0u32;
|
||
for o in obs {
|
||
let key = (
|
||
(o.lat * 1000.0).round() as i32,
|
||
(o.lon * 1000.0).round() as i32,
|
||
);
|
||
// merged's keys are whatever extract_grid produced; use direct
|
||
// lat/lon match by scanning. This is O(N*M) worst-case but N is
|
||
// ~92k and obs mirrors the same point list, so in practice each
|
||
// observation matches exactly one cell by int-key.
|
||
if let Some(cell) = grid.get_mut(&key) {
|
||
if o.max_reflectivity_dbz > 0.0 {
|
||
with_echo += 1;
|
||
}
|
||
cell.insert(
|
||
"nexrad_max_reflectivity_dbz".into(),
|
||
o.max_reflectivity_dbz as f32,
|
||
);
|
||
}
|
||
}
|
||
with_echo
|
||
}
|
||
|
||
fn apply_commercial(grid: &mut PointGrid, lookup: &[LinkLookupEntry]) -> u32 {
|
||
let mut boosted = 0u32;
|
||
for (key, cell) in grid.iter_mut() {
|
||
let (lat, lon) = decoder::key_to_latlon(*key);
|
||
if let Some(d) = commercial::degradation_at((lat, lon), lookup) {
|
||
cell.insert("commercial_degradation_db".into(), d.degradation_db as f32);
|
||
cell.insert("commercial_baseline_dbm".into(), d.baseline_dbm as f32);
|
||
cell.insert("commercial_current_dbm".into(), d.current_dbm as f32);
|
||
cell.insert("commercial_n_links".into(), d.n_links as f32);
|
||
boosted += 1;
|
||
}
|
||
}
|
||
boosted
|
||
}
|
||
|
||
/// Turn one enriched cell into a ProfilesFile entry. The msgpack
|
||
/// `profile` map keeps string keys; the Elixir reader atomizes the
|
||
/// whitelisted subset (surface_*, hpbl_m, pwat_mm, duct metrics, …).
|
||
fn cell_to_profile_entry(lat: f64, lon: f64, cell: &CellValues) -> CellEntry {
|
||
use profiles_file::value_map;
|
||
use rmpv::Value as V;
|
||
|
||
let mut kv: Vec<(&'static str, V)> = Vec::with_capacity(16);
|
||
|
||
if let Some(&v) = cell.get("TMP:2 m above ground") {
|
||
kv.push(("surface_temp_c", V::F64((v as f64) - 273.15)));
|
||
}
|
||
if let Some(&v) = cell.get("DPT:2 m above ground") {
|
||
kv.push(("surface_dewpoint_c", V::F64((v as f64) - 273.15)));
|
||
}
|
||
if let Some(&v) = cell.get("PRES:surface") {
|
||
kv.push(("surface_pressure_mb", V::F64((v as f64) / 100.0)));
|
||
}
|
||
if let Some(&v) = cell.get("HPBL:surface") {
|
||
kv.push(("hpbl_m", V::F64(v as f64)));
|
||
}
|
||
if let Some(&v) = cell.get("PWAT:entire atmosphere (considered as a single layer)") {
|
||
kv.push(("pwat_mm", V::F64(v as f64)));
|
||
}
|
||
// Wind / cloud / precip are read by Elixir's per-QSO scorer
|
||
// (`Propagation.factors_for` in propagation.ex) when a user clicks a
|
||
// grid cell on /map. Without these, wind/sky/rain factors silently
|
||
// score 0 for Rust-produced ProfilesFile cells.
|
||
if let Some(&u) = cell.get("UGRD:10 m above ground") {
|
||
kv.push(("wind_u", V::F64(u as f64)));
|
||
}
|
||
if let Some(&vv) = cell.get("VGRD:10 m above ground") {
|
||
kv.push(("wind_v", V::F64(vv as f64)));
|
||
}
|
||
if let Some(&c) = cell.get("TCDC:entire atmosphere") {
|
||
kv.push(("cloud_cover_pct", V::F64(c as f64)));
|
||
}
|
||
if let Some(&p) = cell.get("APCP:surface") {
|
||
kv.push(("precip_mm", V::F64(p as f64)));
|
||
}
|
||
if let Some(&v) = cell.get("native_min_gradient") {
|
||
kv.push(("native_min_gradient", V::F64(v as f64)));
|
||
}
|
||
if let Some(&v) = cell.get("best_duct_freq_ghz") {
|
||
kv.push(("best_duct_freq_ghz", V::F64(v as f64)));
|
||
}
|
||
if let Some(&v) = cell.get("max_duct_thickness_m") {
|
||
kv.push(("max_duct_thickness_m", V::F64(v as f64)));
|
||
}
|
||
if let Some(&v) = cell.get("duct_count") {
|
||
kv.push(("duct_count", V::Integer((v as i64).into())));
|
||
}
|
||
if let Some(&v) = cell.get("nexrad_max_reflectivity_dbz") {
|
||
kv.push(("nexrad_max_reflectivity_dbz", V::F64(v as f64)));
|
||
}
|
||
if let Some(&d) = cell.get("commercial_degradation_db") {
|
||
let links = cell.get("commercial_n_links").copied().unwrap_or(0.0);
|
||
let baseline = cell.get("commercial_baseline_dbm").copied().unwrap_or(0.0);
|
||
let current = cell.get("commercial_current_dbm").copied().unwrap_or(0.0);
|
||
kv.push((
|
||
"commercial_link_degradation",
|
||
value_map([
|
||
("degradation_db", V::F64(d as f64)),
|
||
("baseline_dbm", V::F64(baseline as f64)),
|
||
("current_dbm", V::F64(current as f64)),
|
||
("n_links", V::Integer((links as i64).into())),
|
||
]),
|
||
));
|
||
}
|
||
|
||
// Pressure-level profile list for SoundingParams.derive on the
|
||
// Elixir side. Mirrors the fields derive expects.
|
||
// Also collect typed `Level` structs so we can pre-compute
|
||
// surface_refractivity and min_refractivity_gradient — the Elixir
|
||
// Skew-T page needs these but SoundingParams.derive will produce nil
|
||
// when any level is missing dewpoint. Pre-computing in Rust (which
|
||
// filters the same way) gives the Elixir side a guaranteed fallback.
|
||
let mut typed_levels: Vec<Level> = Vec::with_capacity(fetcher::grid_level_keys().len());
|
||
let levels: Vec<V> = fetcher::grid_level_keys()
|
||
.iter()
|
||
.filter_map(|k| {
|
||
let t = cell.get(k.tmp.as_str()).copied()?;
|
||
let d = cell.get(k.dpt.as_str()).copied();
|
||
let h = cell.get(k.hgt.as_str()).copied()?;
|
||
let tmpc = (t as f64) - 273.15;
|
||
let dwpc = d.map(|dv| (dv as f64) - 273.15);
|
||
typed_levels.push(Level {
|
||
pres_mb: k.pres_mb,
|
||
hght_m: h as f64,
|
||
tmpc,
|
||
dwpc,
|
||
});
|
||
let mut pairs: Vec<(V, V)> = Vec::with_capacity(4);
|
||
pairs.push((V::String("pres_mb".into()), V::F64(k.pres_mb)));
|
||
pairs.push((V::String("hght_m".into()), V::F64(h as f64)));
|
||
pairs.push((V::String("tmpc".into()), V::F64(tmpc)));
|
||
if let Some(dv) = dwpc {
|
||
pairs.push((V::String("dwpc".into()), V::F64(dv)));
|
||
}
|
||
Some(V::Map(pairs))
|
||
})
|
||
.collect();
|
||
if !levels.is_empty() {
|
||
kv.push(("profile", V::Array(levels)));
|
||
}
|
||
|
||
// Pre-compute refractivity scalars so Elixir consumers (Skew-T,
|
||
// contact detail) have a guaranteed fallback when the pressure-level
|
||
// profile is missing dewpoint on some levels.
|
||
if let Some(sr) = sounding_params::surface_refractivity(&typed_levels) {
|
||
kv.push(("surface_refractivity", V::F64(sr)));
|
||
}
|
||
if let Some(mg) = sounding_params::min_refractivity_gradient(typed_levels) {
|
||
kv.push(("min_refractivity_gradient", V::F64(mg)));
|
||
}
|
||
|
||
let profile = value_map(kv);
|
||
|
||
CellEntry { lat, lon, profile }
|
||
}
|