feat(rust/hrdps): prop-grid-rs HRDPS branch

Stages 4-5 of the HRDPS plan. Lands the Rust forecast-hour pipeline
for source='hrdps' grid_tasks rows.

What ships:

- `hrdps_fetcher.rs` — sibling of fetcher.rs for ECCC's MSC Datamart
  shape: date-prefixed URLs, one variable per GRIB2 file, concurrent
  fetch + byte-concat into a single multi-record blob. wgrib2 reads
  multi-record files natively so the downstream decoder path is
  unchanged. Mirrors lib/microwaveprop/weather/hrdps_client.ex.
- `grid::hrdps_only_points()` + `grid::hrdps_grid_spec()` — Canadian
  cells (49-60°N, -141 to -52°W) minus the HRRR overlap. Disjoint
  from conus_points by construction.
- `db::TaskSource` enum + `source` column on ClaimedTask. claim_next
  surfaces source for forecast lane dispatch; claim_next_analysis
  filters to source='hrrr' only (HRDPS analysis path is not yet
  ported — would need native-duct + NEXRAD + commercial merges that
  Rust doesn't have for HRDPS).
- `pipeline::run_chain_step_hrdps` — fetch combined blob, decode with
  HRDPS grid spec, drop HRRR-overlap cells via a HashSet of
  hrdps_only_points, back-fill DPT from DEPR (HRDPS publishes DEPR
  as primary), score every band, write `<band>/<iso>.hrdps.prop`.
  Skips per-valid_time profile + scalar artifacts since those would
  clobber HRRR's at the same valid_time — surfacing HRDPS on /weather
  is a separate stage.
- worker.rs dispatches forecast lane on (kind, source) — Hrrr →
  run_chain_step, Hrdps → run_chain_step_hrdps. Analysis stays
  HRRR-only. Worker constructs an HrdpsClient alongside HrrrClient.
- `scores_file::write_atomic_hrdps` writes to `.hrdps.prop` with the
  Canadian grid spec via the new encode_with_spec. Coexists with the
  HRRR `.prop` file for the same (band, valid_time).

163 unit tests pass; backfill_dpt_from_depr handles surface and
pressure-level levels with no double-fill when DPT exists.
This commit is contained in:
Graham McIntire 2026-04-29 17:29:16 -05:00
parent fe1c10ce0c
commit 46c2d84e5c
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
7 changed files with 811 additions and 21 deletions

View file

@ -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<dyn std::error::Error>> {
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<dyn std::error::Error>> {
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<HrrrClient>,
hrdps_client: Arc<HrdpsClient>,
scores_dir: PathBuf,
shutdown: Arc<AtomicBool>,
) {
@ -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,

View file

@ -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<Utc>,
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<Option<ClaimedTask>, 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<Option<ClaimedTask>, 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<Option<ClaimedTask>, DbError> {
valid_time,
attempt: attempt + 1,
kind,
source,
}))
}
@ -145,13 +181,18 @@ pub async fn claim_next(pool: &PgPool) -> Result<Option<ClaimedTask>, 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<Option<ClaimedTask>, 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<Option<ClaimedTask>, 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<Option<ClaimedTask>, 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();

View file

@ -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);
}
}

View file

@ -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<Utc>,
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<String>) -> Result<Self, HrdpsFetchError> {
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, HrdpsFetchError> {
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<Utc>,
forecast_hour: u8,
) -> Result<Vec<u8>, HrdpsFetchError> {
let manifest = variable_manifest();
let mut futs: JoinSet<Result<(usize, Vec<u8>), 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<u8>)> = 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<Vec<u8>, 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<Utc>) -> DateTime<Utc> {
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);
}
}

View file

@ -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;

View file

@ -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 `<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);
}
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<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();
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<str>, 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<str> = format!("DPT:{}", depr_key).into();
if cell.contains_key(&dpt_key) {
return None;
}
let tmp_lookup: std::sync::Arc<str> = 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

View file

@ -34,9 +34,11 @@ pub const NO_DATA: u8 = 255;
static UNIQUE: AtomicU64 = AtomicU64::new(1);
/// Absolute path `<base>/<band_mhz>/<iso>.prop`. New writes use this;
/// readers additionally accept the legacy `.ntms` extension while
/// previously-written files drain through the normal retention window.
/// Absolute path `<base>/<band_mhz>/<iso>.prop`. New HRRR writes use
/// this; HRDPS writes go to `<base>/<band_mhz>/<iso>.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<Utc>) -> 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<Utc>) -> Pa
.join(format!("{iso}.prop"))
}
/// HRDPS-source path `<base>/<band_mhz>/<iso>.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<Utc>) -> 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<Utc>) -> PathBuf {
@ -61,14 +74,29 @@ pub enum WriteError {
GridTooLarge { rows: usize, cols: usize },
}
/// Encode the full on-disk representation into a `Vec<u8>`. Separated
/// from the file I/O so tests can round-trip without a tempdir.
/// Encode the full on-disk representation into a `Vec<u8>` 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<Utc>,
scores: &[ScorePoint],
) -> Result<Vec<u8>, 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<Utc>,
scores: &[ScorePoint],
spec: grid::GridSpec,
) -> Result<Vec<u8>, WriteError> {
let n_rows: u16 = spec
.lat_count
.try_into()
@ -111,16 +139,35 @@ pub fn encode(
Ok(out)
}
/// Atomic write to `<base>/<band_mhz>/<iso>.prop`. Creates the band
/// subdirectory on demand.
/// Atomic write to `<base>/<band_mhz>/<iso>.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<Utc>,
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 `<base>/<band_mhz>/<iso>.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<Utc>,
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<u8>) -> 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,