prop/rust/prop_grid_rs/src/hrdps_fetcher.rs
Graham McIntire 63f25a9612
Some checks failed
Build prop-grid-rs / Test, build, push (push) Successful in 5m54s
Build and Push / Build and Push Docker Image (push) Failing after 4m44s
perf(grid-rs): dense grid, fused scoring pass, and columnar .pgrid profiles
Reworks the post-fetch half of the propagation pipeline. Fetch and GRIB2
decode were already cheap — measured against a live HRRR cycle, all 39
pressure messages decode via `wgrib2 -lola` in 0.29 s and the 31 MB
byte-range fetch takes ~3 s — so nothing here touches the decoder. All the
cost was downstream.

Also fixes a broken NOTIFY that made every chain step run up to 5 times.

pg_notify
  `NOTIFY propagation_ready, $1` is a Postgres syntax error: NOTIFY is a
  utility statement whose payload must be a literal, so a bind raises
  42601. It shared a transaction with the `status='done'` UPDATE, so every
  successful step rolled back, stayed 'running', and was requeued by
  reclaim_stale_running up to @max_reclaim_attempts times. Elixir's
  NotifyListener never fired either, so ScoreCache warm and the
  "propagation:updated" fan-out were dead.

FieldGrid
  A decoded grid was HashMap<(i32,i32), HashMap<Arc<str>, f32>> — a dense
  rectangular grid stored as ~95k nested hash maps, costing ~4.6M inserts
  on decode, ~3.7M on merge and ~14M lookups across three derivation
  passes. wgrib2 -lola already emits one dense row-major f32 block per
  message, so keep it: dense per-message planes, names hashed once per
  grid into plane ids, NaN as the missing sentinel. This is what forced
  PROP_GRID_RS_PARALLELISM=1 under a 3Gi limit.

Fused pass
  Three 95k-cell derivation passes plus 23 band-major scoring passes over
  a staged Vec<(f64,f64,Conditions,BandInvariants)> (~19MB re-streamed 23
  times) collapse into one pass: levels extracted once per cell, all 23
  bands scored while the cell is hot, scores accumulated cell-major so
  rayon chunks own disjoint slices. Scores land straight in the dense
  score-file body — no ScorePoint scatter.

.pgrid
  The profile artifact was an rmpv tree plus gzip -9, written 30x an hour,
  and ProfilesFile.read_point/3 gunzipped and unpacked the entire 95k-cell
  file to return one cell on every map click and Skew-T load. Replaced
  with a dense cell-major f32 record array carrying a self-describing
  field table. Elixir reads it via :file.pread; .mp.gz and .etf.gz remain
  readable so files written before this drain out of the 48h window.

  Measured on a full CONUS grid (95,073 cells x 48 planes x 23 bands):
    derive + score + build artifacts   0.022 s
    profile write   3.957 s -> 0.006 s (22.0 MB -> 22.4 MB on disk)
    single-cell read   whole-file decode -> 0.5 us
    23 score files     0.003 s

Also
  - hrrr_points: batched UNNEST upsert replacing one awaited INSERT per
    point. Keeps ON CONFLICT DO UPDATE — the PSKR sampler's two-pass loop
    depends on it.
  - fetcher: real semaphore capping in-flight ranges at
    MAX_PARALLEL_RANGES, which the comment claimed but the code did not do
    (it spawned all 27 while the connection pool was sized for 8).
  - metrics: per-stage histogram. Only chain-step and decode durations
    were instrumented, which is why the write cost stayed invisible.
  - profiles_file: parse_valid_time anchors on the known extension set, so
    sibling-suffixed names like <iso>.hrdps.prop no longer parse as
    <iso>.hrdps and vanish from prune and list operations.
  - PROP_GRID_RS_PARALLELISM 1 -> 3. Memory limit held at 3Gi until RSS is
    observed at the new parallelism.
  - cargo fmt over the crate; worker.rs, hrdps_fetcher.rs and nexrad.rs
    were already unformatted at HEAD and the pre-commit hook gates on it.

HRDPS still runs at 0.5 degrees. wgrib2 -lola scales linearly in output
points on rotated lat/lon (12.5 s wall, 202 s CPU for one message at
0.125 degrees) because it has no inverse projection for those grids; a raw
native dump is 0.32 s. The fix is decode-once plus a closed-form
rotated-pole index, left for a follow-up.
2026-08-01 08:23:36 -05:00

332 lines
12 KiB
Rust
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

//! 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.
match joined {
Ok(inner) => results.push(inner?),
Err(e) => {
tracing::error!(
"HRDPS fetch task panicked: {} ({} remaining variable fetches dropped)",
e,
futs.len()
);
// Drain and collect whatever succeeded.
while let Some(rest) = futs.join_next().await {
if let Ok(Ok((idx, buf))) = rest {
results.push((idx, buf));
}
}
break;
}
}
}
// 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);
}
}