diff --git a/rust/prop_grid_rs/src/bin/worker.rs b/rust/prop_grid_rs/src/bin/worker.rs index 9dfad794..3dfb69a5 100644 --- a/rust/prop_grid_rs/src/bin/worker.rs +++ b/rust/prop_grid_rs/src/bin/worker.rs @@ -16,7 +16,7 @@ use std::sync::atomic::{AtomicBool, Ordering}; use std::sync::Arc; use std::time::{Duration, Instant}; -use prop_grid_rs::{db, fetcher::HrrrClient, metrics, pipeline, telemetry}; +use prop_grid_rs::{db, fetcher::HrrrClient, hrdps_fetcher::HrdpsClient, metrics, pipeline, telemetry}; use tokio::task::JoinSet; use tracing::{error, info, warn}; @@ -63,6 +63,12 @@ async fn main() -> Result<(), Box> { let pool = db::connect(&database_url, max_conn).await?; let client = Arc::new(HrrrClient::new_with_fallback(&hrrr_base, &hrrr_fallback)?); + // HRDPS shares the same worker pool. The MSC datamart base rarely + // changes; allow override only via env so production doesn't need a + // redeploy if ECCC moves the public mirror. + let hrdps_base = std::env::var("HRDPS_DATAMART_BASE") + .unwrap_or_else(|_| prop_grid_rs::hrdps_fetcher::DATAMART_BASE_DEFAULT.to_string()); + let hrdps_client = Arc::new(HrdpsClient::new(&hrdps_base)?); let shutdown = Arc::new(AtomicBool::new(false)); spawn_signal_handler(shutdown.clone()); @@ -89,10 +95,11 @@ async fn main() -> Result<(), Box> { for worker_id in 0..parallelism { let pool = pool.clone(); let client = client.clone(); + let hrdps_client = hrdps_client.clone(); let scores_dir = scores_dir.clone(); let shutdown = shutdown.clone(); set.spawn(async move { - worker_loop(worker_id, pool, client, scores_dir, shutdown).await; + worker_loop(worker_id, pool, client, hrdps_client, scores_dir, shutdown).await; }); } @@ -113,6 +120,7 @@ async fn worker_loop( worker_id: usize, pool: sqlx::PgPool, client: Arc, + hrdps_client: Arc, scores_dir: PathBuf, shutdown: Arc, ) { @@ -133,8 +141,8 @@ async fn worker_loop( let task_id = task.id; let _guard = metrics::InFlightGuard::new(); let started = Instant::now(); - let result = match task.kind { - db::TaskKind::Forecast => pipeline::run_chain_step( + let result = match (task.kind.clone(), task.source) { + (db::TaskKind::Forecast, db::TaskSource::Hrrr) => pipeline::run_chain_step( &client, &scores_dir, &pipeline::ChainStepInput { @@ -144,7 +152,19 @@ async fn worker_loop( ) .await .map(ChainOutcome::Forecast), - db::TaskKind::Analysis => pipeline::run_analysis_step( + (db::TaskKind::Forecast, db::TaskSource::Hrdps) => { + pipeline::run_chain_step_hrdps( + &hrdps_client, + &scores_dir, + &pipeline::ChainStepInput { + run_time: task.run_time, + forecast_hour: task.forecast_hour as u8, + }, + ) + .await + .map(ChainOutcome::Forecast) + } + (db::TaskKind::Analysis, _) => pipeline::run_analysis_step( &client, &pool, &scores_dir, @@ -165,6 +185,7 @@ async fn worker_loop( worker_id, task_id = %task.id, kind = "forecast", + source = task.source.as_str(), run_time = %task.run_time, forecast_hour = task.forecast_hour, files = stats.score_files_written, diff --git a/rust/prop_grid_rs/src/db.rs b/rust/prop_grid_rs/src/db.rs index e55753e9..d23d953a 100644 --- a/rust/prop_grid_rs/src/db.rs +++ b/rust/prop_grid_rs/src/db.rs @@ -39,6 +39,38 @@ impl TaskKind { } } +/// NWP source the task should fetch from. Set by the Elixir seeder +/// (`GridTaskEnqueuer.seed_with_analysis/2`); the worker dispatches the +/// fetch + decode + score pipeline on this enum. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum TaskSource { + /// HRRR — CONUS coverage, hourly cycles, NOAA AWS S3. + Hrrr, + /// HRDPS — Canadian extent (49-60°N, -141 to -52°W minus the HRRR + /// overlap), 4×/day at 00/06/12/18Z, MSC Datamart. + Hrdps, +} + +impl TaskSource { + fn from_str(s: &str) -> Self { + match s { + "hrdps" => TaskSource::Hrdps, + // Default to Hrrr for the legacy rows that predated the + // source column (and for any source string we don't + // recognize — we'd rather keep doing the safe HRRR path + // than refuse to claim). + _ => TaskSource::Hrrr, + } + } + + pub fn as_str(self) -> &'static str { + match self { + TaskSource::Hrrr => "hrrr", + TaskSource::Hrdps => "hrdps", + } + } +} + #[derive(Debug, Clone)] pub struct ClaimedTask { pub id: Uuid, @@ -47,6 +79,7 @@ pub struct ClaimedTask { pub valid_time: DateTime, pub attempt: i32, pub kind: TaskKind, + pub source: TaskSource, } #[derive(Debug, thiserror::Error)] @@ -82,7 +115,7 @@ pub async fn claim_next(pool: &PgPool) -> Result, DbError> { // reseeds the new run and the prune cron clears rotted rows. let row = sqlx::query( r#" - SELECT id, run_time, forecast_hour, valid_time, attempt, kind + SELECT id, run_time, forecast_hour, valid_time, attempt, kind, source FROM grid_tasks WHERE status = 'queued' AND forecast_hour > 0 AND kind = 'forecast' ORDER BY run_time DESC, forecast_hour ASC @@ -112,6 +145,8 @@ pub async fn claim_next(pool: &PgPool) -> Result, DbError> { let attempt: i32 = row.try_get("attempt")?; let kind_str: String = row.try_get("kind")?; let kind = TaskKind::from_str(&kind_str); + let source_str: String = row.try_get("source")?; + let source = TaskSource::from_str(&source_str); sqlx::query( r#" @@ -136,6 +171,7 @@ pub async fn claim_next(pool: &PgPool) -> Result, DbError> { valid_time, attempt: attempt + 1, kind, + source, })) } @@ -145,13 +181,18 @@ pub async fn claim_next(pool: &PgPool) -> Result, DbError> { /// learned to do in Phase 3 Stream A. Kept on a separate lane from /// `claim_next` so operators can gate the analysis path independently /// of the forecast path during rollout. +/// +/// Restricted to `source = 'hrrr'` for now — the HRDPS analysis path +/// (native-level duct + NEXRAD + commercial merges) hasn't been ported +/// yet, so HRDPS analysis rows would fail mid-pipeline. HRDPS analysis +/// stays Elixir-side / out of scope until that gap is closed. pub async fn claim_next_analysis(pool: &PgPool) -> Result, DbError> { let mut tx = pool.begin().await?; let row = sqlx::query( r#" - SELECT id, run_time, forecast_hour, valid_time, attempt, kind + SELECT id, run_time, forecast_hour, valid_time, attempt, kind, source FROM grid_tasks - WHERE status = 'queued' AND kind = 'analysis' + WHERE status = 'queued' AND kind = 'analysis' AND source = 'hrrr' ORDER BY run_time DESC, forecast_hour ASC FOR UPDATE SKIP LOCKED LIMIT 1 @@ -175,6 +216,8 @@ pub async fn claim_next_analysis(pool: &PgPool) -> Result, D let attempt: i32 = row.try_get("attempt")?; let kind_str: String = row.try_get("kind")?; let kind = TaskKind::from_str(&kind_str); + let source_str: String = row.try_get("source")?; + let source = TaskSource::from_str(&source_str); sqlx::query( r#" @@ -199,6 +242,7 @@ pub async fn claim_next_analysis(pool: &PgPool) -> Result, D valid_time, attempt: attempt + 1, kind, + source, })) } @@ -416,6 +460,8 @@ mod tests { assert!(claimed.forecast_hour > 0); assert_eq!(claimed.attempt, 1); assert_eq!(claimed.kind, TaskKind::Forecast); + // Default-DB rows land with source='hrrr' per the migration default. + assert_eq!(claimed.source, TaskSource::Hrrr); complete(&pool, &claimed).await.unwrap(); diff --git a/rust/prop_grid_rs/src/grid.rs b/rust/prop_grid_rs/src/grid.rs index 65d6483a..c1a70591 100644 --- a/rust/prop_grid_rs/src/grid.rs +++ b/rust/prop_grid_rs/src/grid.rs @@ -1,5 +1,13 @@ -//! CONUS grid at 0.125° resolution. 1:1 port of -//! `lib/microwaveprop/propagation/grid.ex`. +//! Grid definitions at 0.125° resolution. 1:1 port of +//! `lib/microwaveprop/propagation/grid.ex`. Two regions, disjoint by +//! construction: +//! +//! * `conus_points()` — HRRR coverage (lat 25-50, lon -125 to -66). +//! * `hrdps_only_points()` — Canadian extent covered by HRDPS but not +//! HRRR (lat 49-60 minus the HRRR overlap, lon -141 to -52). +//! +//! Coverage stops at 60°N for v1 (SRTM elevation cuts off there; Arctic +//! CDEM coverage is deferred to a follow-up plan). pub const LAT_MIN: f64 = 25.0; pub const LAT_MAX: f64 = 50.0; @@ -7,6 +15,11 @@ pub const LON_MIN: f64 = -125.0; pub const LON_MAX: f64 = -66.0; pub const STEP: f64 = 0.125; +pub const HRDPS_LAT_MIN: f64 = 49.0; +pub const HRDPS_LAT_MAX: f64 = 60.0; +pub const HRDPS_LON_MIN: f64 = -141.0; +pub const HRDPS_LON_MAX: f64 = -52.0; + #[derive(Debug, Clone, Copy, PartialEq)] pub struct GridSpec { pub lon_start: f64, @@ -43,6 +56,47 @@ pub fn conus_points() -> Vec<(f64, f64)> { out } +/// HRDPS grid spec covering the full Canadian bbox. The `hrdps_only_points()` +/// list is a subset (excludes the HRRR overlap), but the underlying wgrib2 +/// extraction needs the full bbox so the score-file's lat_start/lon_start +/// align with the cells the worker actually wrote. +pub fn hrdps_grid_spec() -> GridSpec { + let lon_count = ((HRDPS_LON_MAX - HRDPS_LON_MIN) / STEP).round() as usize + 1; + let lat_count = ((HRDPS_LAT_MAX - HRDPS_LAT_MIN) / STEP).round() as usize + 1; + GridSpec { + lon_start: HRDPS_LON_MIN, + lon_count, + lon_step: STEP, + lat_start: HRDPS_LAT_MIN, + lat_count, + lat_step: STEP, + } +} + +/// Canadian-only grid points: cells inside the HRDPS bbox (49-60°N, +/// -141 to -52°W) but outside HRRR's CONUS bbox. Disjoint from +/// `conus_points()` by construction so the two grids never double-write +/// the same `(lat, lon)`. +pub fn hrdps_only_points() -> Vec<(f64, f64)> { + let spec = hrdps_grid_spec(); + let mut out = Vec::with_capacity(spec.lat_count * spec.lon_count); + for j in 0..spec.lat_count { + let lat = round3(spec.lat_start + j as f64 * spec.lat_step); + for i in 0..spec.lon_count { + let lon = round3(spec.lon_start + i as f64 * spec.lon_step); + if !in_conus_bbox(lat, lon) { + out.push((lat, lon)); + } + } + } + out +} + +#[inline] +fn in_conus_bbox(lat: f64, lon: f64) -> bool { + lat >= LAT_MIN && lat <= LAT_MAX && lon >= LON_MIN && lon <= LON_MAX +} + /// Matches Elixir's `Float.round(x, 3)` — banker's rounding isn't used, /// half-away-from-zero is. For grid coords this matches BEAM's behavior /// for the values we care about (no exact half cases). @@ -79,4 +133,46 @@ mod tests { assert_eq!(points[0], (25.0, -125.0)); assert_eq!(points[points.len() - 1], (50.0, -66.0)); } + + #[test] + fn hrdps_grid_spec_matches_elixir_bbox() { + let spec = hrdps_grid_spec(); + // Elixir Grid: hrdps_lat_min=49, hrdps_lat_max=60, hrdps_lon_min=-141, hrdps_lon_max=-52 + // lon_count = (-52 - -141) / 0.125 + 1 = 713 + // lat_count = (60 - 49) / 0.125 + 1 = 89 + assert_eq!(spec.lon_count, 713); + assert_eq!(spec.lat_count, 89); + assert_eq!(spec.lon_start, -141.0); + assert_eq!(spec.lat_start, 49.0); + } + + #[test] + fn hrdps_only_points_are_disjoint_from_conus() { + let conus: std::collections::HashSet<(u64, u64)> = conus_points() + .into_iter() + .map(|(la, lo)| (la.to_bits(), lo.to_bits())) + .collect(); + + for (la, lo) in hrdps_only_points() { + assert!( + !conus.contains(&(la.to_bits(), lo.to_bits())), + "hrdps point {la},{lo} overlaps conus" + ); + } + } + + #[test] + fn hrdps_only_points_all_within_bbox() { + for (lat, lon) in hrdps_only_points() { + assert!((49.0..=60.0).contains(&lat)); + assert!((-141.0..=-52.0).contains(&lon)); + } + } + + #[test] + fn hrdps_only_includes_cells_above_50n() { + // At least one cell strictly north of HRRR's lat_max=50. + let any_above_50 = hrdps_only_points().iter().any(|(la, _)| *la > 50.0); + assert!(any_above_50); + } } diff --git a/rust/prop_grid_rs/src/hrdps_fetcher.rs b/rust/prop_grid_rs/src/hrdps_fetcher.rs new file mode 100644 index 00000000..069f73ee --- /dev/null +++ b/rust/prop_grid_rs/src/hrdps_fetcher.rs @@ -0,0 +1,312 @@ +//! HRDPS (Canadian) fetcher. Sibling of `fetcher.rs` (HRRR) but with +//! ECCC's MSC Datamart shape: +//! +//! * Date-prefixed URL structure +//! (`https://dd.weather.gc.ca/{YYYYMMDD}/WXO-DD/model_hrdps/...`) +//! * **One variable per GRIB2 file** (~3.5 MB each, ~14 files total). +//! No idx / no byte-range partials. +//! * Concurrent per-variable fetch + byte-concatenate into a single +//! multi-record blob — wgrib2 reads multi-record files natively, so +//! the downstream `decoder::extract_grid` path is unchanged from HRRR. +//! * Rotated lat/lon projection handled transparently by wgrib2. +//! +//! 1:1 port of `lib/microwaveprop/weather/hrdps_client.ex`. + +use std::time::Duration; + +use chrono::{DateTime, Timelike, Utc}; +use reqwest::Client; +use tokio::task::JoinSet; + +pub const DATAMART_BASE_DEFAULT: &str = "https://dd.weather.gc.ca"; +pub const REQUEST_TIMEOUT: Duration = Duration::from_secs(120); +pub const MAX_PARALLEL_FETCHES: usize = 4; + +/// Pressure levels matching `lib/microwaveprop/weather/hrdps_client.ex`'s +/// `@grid_pressure_levels`. Mirrors HRRR's grid set (1000-700 mb) so +/// SoundingParams.derive sees a comparable refractivity-gradient profile. +pub const HRDPS_GRID_PRESSURE_LEVELS: &[u16] = &[1000, 950, 900, 850, 800, 750, 700]; + +/// Surface variables we fetch per cycle/forecast hour. Each entry is +/// (internal_atom, msc_var, msc_level) — the MSC level slug becomes part +/// of the filename. Mirrors `HrdpsClient.@surface_vars` in Elixir. +pub const HRDPS_SURFACE_VARS: &[(&str, &str, &str)] = &[ + ("tmp_2m", "TMP", "AGL-2m"), + ("depr_2m", "DEPR", "AGL-2m"), + ("dpt_2m", "DPT", "AGL-2m"), + ("pres_sfc", "PRES", "Sfc"), + ("hpbl_sfc", "HPBL", "Sfc"), + ("ugrd_10m", "UGRD", "AGL-10m"), + ("vgrd_10m", "VGRD", "AGL-10m"), + ("tcdc_sfc", "TCDC", "Sfc"), +]; + +#[derive(Debug, thiserror::Error)] +pub enum HrdpsFetchError { + #[error("http: {0}")] + Http(#[from] reqwest::Error), + #[error("HRDPS file fetch HTTP {status} for {url}")] + Status { status: u16, url: String }, + #[error("HRDPS fetch returned empty body for {url}")] + EmptyBody { url: String }, +} + +/// Build the MSC Datamart URL for one HRDPS variable file. ECCC reorganized +/// to a date-prefixed structure as of 2026 — the old `/model_hrdps/...` +/// flat root 404s. Filename pattern verified against the live datamart on +/// 2026-04-29. +pub fn hrdps_url( + base: &str, + cycle: DateTime, + forecast_hour: u8, + msc_var: &str, + msc_level: &str, +) -> String { + let date = cycle.date_naive(); + let date_str = date.format("%Y%m%d").to_string(); + let hour_str = format!("{:02}", cycle.hour()); + let fff_str = format!("{:03}", forecast_hour); + let trimmed = base.trim_end_matches('/'); + + format!( + "{trimmed}/{date_str}/WXO-DD/model_hrdps/continental/2.5km/{hour_str}/{fff_str}/\ + {date_str}T{hour_str}Z_MSC_HRDPS_{msc_var}_{msc_level}_RLatLon0.0225_PT{fff_str}H.grib2" + ) +} + +/// Build the variable manifest the fetcher pulls per cycle/forecast hour. +/// Returns `[(msc_var, msc_level)]` in stable order so cargo-test fixtures +/// match production output byte-for-byte. Surface variables come first, +/// then pressure-level (TMP / DEPR / HGT × levels). +pub fn variable_manifest() -> Vec<(&'static str, String)> { + let mut out: Vec<(&'static str, String)> = HRDPS_SURFACE_VARS + .iter() + .map(|(_, var, level)| (*var, level.to_string())) + .collect(); + + for &level in HRDPS_GRID_PRESSURE_LEVELS { + let level_slug = format!("ISBL_{:04}", level); + out.push(("TMP", level_slug.clone())); + out.push(("DEPR", level_slug.clone())); + out.push(("HGT", level_slug)); + } + + out +} + +#[derive(Debug, Clone)] +pub struct HrdpsClient { + http: Client, + base: String, +} + +impl HrdpsClient { + pub fn new(base: impl Into) -> Result { + let http = Client::builder() + .timeout(REQUEST_TIMEOUT) + .pool_max_idle_per_host(MAX_PARALLEL_FETCHES) + .build()?; + + Ok(Self { + http, + base: base.into(), + }) + } + + pub fn default_base() -> Result { + Self::new(DATAMART_BASE_DEFAULT) + } + + pub fn base(&self) -> &str { + &self.base + } + + /// Fetch every variable file for `(cycle, forecast_hour)` concurrently, + /// then byte-concatenate into a single multi-record GRIB2 blob. wgrib2 + /// reads multi-record files natively, so the returned blob plugs into + /// `decoder::extract_grid` exactly like an HRRR `wrfsfcf`/`wrfprsf` + /// blob would — same `match` pattern works because HRDPS uses the same + /// NCEP-style level strings (`TMP:2 m above ground`, etc.) at the + /// wgrib2 inventory layer regardless of the MSC filename. + pub async fn fetch_combined_blob( + &self, + cycle: DateTime, + forecast_hour: u8, + ) -> Result, HrdpsFetchError> { + let manifest = variable_manifest(); + let mut futs: JoinSet), HrdpsFetchError>> = JoinSet::new(); + + // Cap concurrent fetches at MAX_PARALLEL_FETCHES (datamart isn't + // aggressively rate-limited but ECCC has explicit guidance to + // limit hammer-rate). Index preserves variable order so the + // concatenated output is deterministic. + let manifest_len = manifest.len(); + for (idx, (msc_var, msc_level)) in manifest.into_iter().enumerate() { + let client = self.http.clone(); + let url = hrdps_url(&self.base, cycle, forecast_hour, msc_var, &msc_level); + futs.spawn(async move { + let bytes = retry_get(&client, &url).await?; + Ok((idx, bytes)) + }); + } + + let mut results: Vec<(usize, Vec)> = Vec::with_capacity(manifest_len); + while let Some(joined) = futs.join_next().await { + // Per CLAUDE.md: never silently swallow a task failure. A + // single missing variable here means a Canadian dead-zone + // we won't notice for hours. + results.push(joined.expect("hrdps fetch task panicked")?); + } + + // Stable ordering by manifest index. wgrib2 doesn't care about + // record order but deterministic concatenation makes byte-level + // golden tests possible. + results.sort_by_key(|(idx, _)| *idx); + + let total: usize = results.iter().map(|(_, b)| b.len()).sum(); + let mut out = Vec::with_capacity(total); + for (_, bytes) in results { + out.extend_from_slice(&bytes); + } + Ok(out) + } +} + +async fn retry_get(client: &Client, url: &str) -> Result, HrdpsFetchError> { + for attempt in 0..5 { + match client.get(url).send().await { + Ok(resp) => { + let status = resp.status().as_u16(); + if status == 200 { + let bytes = resp.bytes().await?; + if bytes.is_empty() { + return Err(HrdpsFetchError::EmptyBody { + url: url.to_string(), + }); + } + return Ok(bytes.to_vec()); + } + if !is_retryable_status(status) { + return Err(HrdpsFetchError::Status { + status, + url: url.to_string(), + }); + } + tracing::warn!(status, attempt, url, "hrdps retryable status"); + } + Err(e) if e.is_connect() || e.is_timeout() => { + tracing::warn!(error = %e, attempt, url, "hrdps transient error"); + } + Err(e) => return Err(HrdpsFetchError::Http(e)), + } + tokio::time::sleep(retry_backoff(attempt)).await; + } + Err(HrdpsFetchError::Status { + status: 0, + url: url.to_string(), + }) +} + +fn is_retryable_status(status: u16) -> bool { + matches!(status, 429 | 500 | 502 | 503 | 504) +} + +fn retry_backoff(attempt: u32) -> Duration { + let base_ms = 1000_u64 << attempt; + let jitter = rand_jitter_ms(); + Duration::from_millis(base_ms + jitter) +} + +// Tiny stdlib-only jitter — avoids pulling in `rand` for one call site. +// Returns 0..=999 ms. +fn rand_jitter_ms() -> u64 { + use std::time::{SystemTime, UNIX_EPOCH}; + let nanos = SystemTime::now() + .duration_since(UNIX_EPOCH) + .map(|d| d.subsec_nanos() as u64) + .unwrap_or(0); + nanos % 1_000 +} + +/// Snap a UTC timestamp to the start of the most recent HRDPS cycle hour +/// (00/06/12/18Z). Always rounds DOWN — never points at a future cycle +/// that may not have published yet. +pub fn nearest_hrdps_cycle(dt: DateTime) -> DateTime { + let cycle_hour = (dt.hour() / 6) * 6; + dt.with_hour(cycle_hour) + .and_then(|d| d.with_minute(0)) + .and_then(|d| d.with_second(0)) + .and_then(|d| d.with_nanosecond(0)) + .unwrap_or(dt) +} + +#[cfg(test)] +mod tests { + use super::*; + use chrono::TimeZone; + + #[test] + fn url_matches_live_datamart_format() { + let cycle = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(); + let url = hrdps_url(DATAMART_BASE_DEFAULT, cycle, 0, "TMP", "AGL-2m"); + assert_eq!( + url, + "https://dd.weather.gc.ca/20260429/WXO-DD/model_hrdps/continental/2.5km/12/000/\ + 20260429T12Z_MSC_HRDPS_TMP_AGL-2m_RLatLon0.0225_PT000H.grib2" + ); + } + + #[test] + fn url_pads_forecast_hour_to_three_digits() { + let cycle = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(); + let url = hrdps_url(DATAMART_BASE_DEFAULT, cycle, 24, "TMP", "AGL-2m"); + assert!(url.contains("/024/")); + assert!(url.contains("PT024H.grib2")); + } + + #[test] + fn url_pads_cycle_hour_to_two_digits() { + let cycle = Utc.with_ymd_and_hms(2026, 4, 29, 6, 0, 0).unwrap(); + let url = hrdps_url(DATAMART_BASE_DEFAULT, cycle, 0, "TMP", "AGL-2m"); + assert!(url.contains("/06/000/")); + assert!(url.contains("T06Z_")); + } + + #[test] + fn pressure_level_url_uses_isbl_padding() { + let cycle = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(); + let url = hrdps_url(DATAMART_BASE_DEFAULT, cycle, 0, "TMP", "ISBL_0850"); + assert!(url.contains("MSC_HRDPS_TMP_ISBL_0850_")); + } + + #[test] + fn manifest_covers_surface_and_pressure_set() { + let manifest = variable_manifest(); + // 8 surface + 7 pressure levels × 3 vars (TMP/DEPR/HGT) = 29 + assert_eq!(manifest.len(), 29); + assert!(manifest.iter().any(|(v, l)| *v == "TMP" && l == "AGL-2m")); + assert!(manifest.iter().any(|(v, l)| *v == "TMP" && l == "ISBL_0850")); + assert!(manifest.iter().any(|(v, l)| *v == "HGT" && l == "ISBL_1000")); + } + + #[test] + fn nearest_cycle_snaps_to_six_hour_boundary() { + let dt = Utc.with_ymd_and_hms(2026, 4, 29, 13, 45, 0).unwrap(); + let cycle = nearest_hrdps_cycle(dt); + assert_eq!(cycle, Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap()); + } + + #[test] + fn nearest_cycle_rounds_down_at_boundary() { + // Just inside the 06Z window; mustn't jump forward to 12Z. + let dt = Utc.with_ymd_and_hms(2026, 4, 29, 11, 59, 0).unwrap(); + let cycle = nearest_hrdps_cycle(dt); + assert_eq!(cycle, Utc.with_ymd_and_hms(2026, 4, 29, 6, 0, 0).unwrap()); + } + + #[test] + fn nearest_cycle_idempotent_on_boundary() { + let dt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(); + assert_eq!(nearest_hrdps_cycle(dt), dt); + } +} diff --git a/rust/prop_grid_rs/src/lib.rs b/rust/prop_grid_rs/src/lib.rs index 691cb531..e9d6b4be 100644 --- a/rust/prop_grid_rs/src/lib.rs +++ b/rust/prop_grid_rs/src/lib.rs @@ -39,6 +39,7 @@ pub mod decoder; pub mod duct; pub mod fetcher; pub mod grid; +pub mod hrdps_fetcher; pub mod hrrr_points; pub mod metrics; pub mod native_duct; diff --git a/rust/prop_grid_rs/src/pipeline.rs b/rust/prop_grid_rs/src/pipeline.rs index 56041f29..f53446c7 100644 --- a/rust/prop_grid_rs/src/pipeline.rs +++ b/rust/prop_grid_rs/src/pipeline.rs @@ -15,7 +15,8 @@ use crate::band_config; use crate::commercial::{self, LinkLookupEntry}; use crate::decoder::{self, CellValues, PointGrid}; use crate::fetcher::{self, HrrrClient, Product}; -use crate::grid::wgrib2_grid_spec; +use crate::grid::{hrdps_grid_spec, 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}; @@ -28,6 +29,8 @@ use crate::weather_scalar_file::{self, ScalarRow}; 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}")] @@ -255,6 +258,170 @@ pub async fn run_chain_step( }) } +/// 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 `//.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 { + 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); + } + + let grid_spec = hrdps_grid_spec(); + // 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):"; + + let pattern_owned = pattern.to_string(); + let blob_for_decode = blob.clone(); + let grid_spec_for_decode = grid_spec; + let mut grid = tokio::task::spawn_blocking(move || { + decoder::extract_grid(&blob_for_decode, &pattern_owned, grid_spec_for_decode) + }) + .await + .expect("blocking join")?; + drop(blob); + + // Restrict to the Canadian-only mask (drop HRRR-overlap cells) + // BEFORE the DEPR fill + scoring loop so the inner work scales + // with ~57k Canadian cells, not 89×713 ≈ 63k Canadian-bbox cells. + let hrdps_keys: std::collections::HashSet<(u64, u64)> = hrdps_only_points() + .into_iter() + .map(|(la, lo)| (la.to_bits(), lo.to_bits())) + .collect(); + + grid.retain(|key, _| { + let (lat, lon) = decoder::key_to_latlon(*key); + hrdps_keys.contains(&(lat.to_bits(), lon.to_bits())) + }); + + let point_count = grid.len() as u32; + + let mut prepared: Vec<(f64, f64, Conditions, scorer::BandInvariants)> = + 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)); + } + } + + let bands = band_config::all_bands(); + let prepared = std::sync::Arc::new(prepared); + let scored: Vec<(u32, Vec)> = { + use rayon::prelude::*; + bands + .par_iter() + .map(|band| { + let mut scores: Vec = 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(); + 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")?; + } + 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, f32)> = cell + .iter() + .filter_map(|(key, _)| { + let depr_key = key.strip_prefix("DEPR:")?; + let tmp_key = format!("TMP:{}", depr_key); + let dpt_key: std::sync::Arc = format!("DPT:{}", depr_key).into(); + if cell.contains_key(&dpt_key) { + return None; + } + let tmp_lookup: std::sync::Arc = tmp_key.into(); + let tmp = cell.get(&tmp_lookup).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 { @@ -470,6 +637,63 @@ mod tests { 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 diff --git a/rust/prop_grid_rs/src/scores_file.rs b/rust/prop_grid_rs/src/scores_file.rs index cd7bc927..aa6c5b86 100644 --- a/rust/prop_grid_rs/src/scores_file.rs +++ b/rust/prop_grid_rs/src/scores_file.rs @@ -34,9 +34,11 @@ pub const NO_DATA: u8 = 255; static UNIQUE: AtomicU64 = AtomicU64::new(1); -/// Absolute path `//.prop`. New writes use this; -/// readers additionally accept the legacy `.ntms` extension while -/// previously-written files drain through the normal retention window. +/// Absolute path `//.prop`. New HRRR writes use +/// this; HRDPS writes go to `//.hrdps.prop` via +/// `path_for_hrdps`. Readers accept both plus the legacy `.ntms` +/// extension while previously-written files drain through the normal +/// retention window. pub fn path_for(base_dir: &Path, band_mhz: u32, valid_time: DateTime) -> PathBuf { let iso = valid_time.format("%Y-%m-%dT%H:%M:%SZ").to_string(); base_dir @@ -44,6 +46,17 @@ pub fn path_for(base_dir: &Path, band_mhz: u32, valid_time: DateTime) -> Pa .join(format!("{iso}.prop")) } +/// HRDPS-source path `//.hrdps.prop`. Sibling of +/// `path_for/3` so HRRR (CONUS) and HRDPS (Canadian) score grids for the +/// same `(band, valid_time)` coexist as separate files. Readers stitch +/// them together at request time. +pub fn path_for_hrdps(base_dir: &Path, band_mhz: u32, valid_time: DateTime) -> PathBuf { + let iso = valid_time.format("%Y-%m-%dT%H:%M:%SZ").to_string(); + base_dir + .join(band_mhz.to_string()) + .join(format!("{iso}.hrdps.prop")) +} + /// Legacy path used by scores written before the `.prop` rename. Kept so /// readers can discover still-on-disk `.ntms` files during rollout. pub fn legacy_path_for(base_dir: &Path, band_mhz: u32, valid_time: DateTime) -> PathBuf { @@ -61,14 +74,29 @@ pub enum WriteError { GridTooLarge { rows: usize, cols: usize }, } -/// Encode the full on-disk representation into a `Vec`. Separated -/// from the file I/O so tests can round-trip without a tempdir. +/// Encode the full on-disk representation into a `Vec` using the +/// CONUS HRRR grid spec. Backwards-compatible wrapper around +/// `encode_with_spec` for HRRR callers; HRDPS callers pass +/// `grid::hrdps_grid_spec()` explicitly. pub fn encode( band_mhz: u32, valid_time: DateTime, scores: &[ScorePoint], ) -> Result, WriteError> { - let spec = grid::wgrib2_grid_spec(); + encode_with_spec(band_mhz, valid_time, scores, grid::wgrib2_grid_spec()) +} + +/// Encode score scrubs into the on-disk format using an explicit +/// `GridSpec`. Used by both HRRR (CONUS spec) and HRDPS (Canadian spec) +/// writers. Cells outside the spec's bbox are silently dropped — for +/// HRDPS that's the expected case for the HRRR-overlap cells we mask +/// out at scoring time anyway. +pub fn encode_with_spec( + band_mhz: u32, + valid_time: DateTime, + scores: &[ScorePoint], + spec: grid::GridSpec, +) -> Result, WriteError> { let n_rows: u16 = spec .lat_count .try_into() @@ -111,16 +139,35 @@ pub fn encode( Ok(out) } -/// Atomic write to `//.prop`. Creates the band -/// subdirectory on demand. +/// Atomic write to `//.prop` (HRRR/CONUS). +/// HRDPS callers use `write_atomic_hrdps` which writes to `.hrdps.prop` +/// with the Canadian grid spec. pub fn write_atomic( base_dir: &Path, band_mhz: u32, valid_time: DateTime, scores: &[ScorePoint], ) -> Result<(), WriteError> { - let bytes = encode(band_mhz, valid_time, scores)?; - let final_path = path_for(base_dir, band_mhz, valid_time); + write_atomic_at( + path_for(base_dir, band_mhz, valid_time), + encode(band_mhz, valid_time, scores)?, + ) +} + +/// Atomic write to `//.hrdps.prop` using the +/// Canadian grid spec. Coexists with the sibling HRRR `.prop` file for +/// the same `(band, valid_time)`. +pub fn write_atomic_hrdps( + base_dir: &Path, + band_mhz: u32, + valid_time: DateTime, + scores: &[ScorePoint], +) -> Result<(), WriteError> { + let bytes = encode_with_spec(band_mhz, valid_time, scores, grid::hrdps_grid_spec())?; + write_atomic_at(path_for_hrdps(base_dir, band_mhz, valid_time), bytes) +} + +fn write_atomic_at(final_path: PathBuf, bytes: Vec) -> Result<(), WriteError> { let parent = final_path .parent() .ok_or_else(|| std::io::Error::other("path has no parent"))?; @@ -152,6 +199,49 @@ pub fn write_atomic( } } +#[cfg(test)] +mod hrdps_tests { + use super::*; + use chrono::TimeZone; + + #[test] + fn hrdps_path_uses_hrdps_prop_extension() { + let base = Path::new("/data/scores"); + let vt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(); + let path = path_for_hrdps(base, 10_000, vt); + assert!(path.to_string_lossy().ends_with("/2026-04-29T12:00:00Z.hrdps.prop")); + } + + #[test] + fn hrdps_encode_uses_canadian_bbox() { + let vt = Utc.with_ymd_and_hms(2026, 4, 29, 12, 0, 0).unwrap(); + let scores = vec![ + ScorePoint { + lat: 53.32, + lon: -60.42, + score: 75, + }, + // Cell outside HRDPS bbox (Hawaii) — should be silently dropped. + ScorePoint { + lat: 21.3, + lon: -157.8, + score: 50, + }, + ]; + + let bytes = + encode_with_spec(10_000, vt, &scores, grid::hrdps_grid_spec()).expect("encodes"); + let decoded = decode(&bytes).expect("decodes"); + + // HRDPS bbox lat_start=49, lon_start=-141 + assert!((decoded.lat_min - 49.0).abs() < 1e-6); + assert!((decoded.lon_min - (-141.0)).abs() < 1e-6); + // 89 lat × 713 lon cells per the Canadian bbox. + assert_eq!(decoded.n_rows, 89); + assert_eq!(decoded.n_cols, 713); + } +} + #[derive(Debug, Clone, Copy)] pub struct ScorePoint { pub lat: f64,