From 2348c48c2634c97dbdf51625551b87289ef680e5 Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Wed, 29 Apr 2026 18:25:40 -0500 Subject: [PATCH] fix(rust/hrdps): replace -lola full-grid extract with -lon point extract MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Production observation 2026-04-29: HRDPS chain steps using decoder::extract_grid (wgrib2 -lola) ran at 1000m CPU for 10+ minutes per task without producing any output. The -lola full-grid interpolation against HRDPS's rotated lat/lon source grid is far slower than against HRRR's Lambert source — wgrib2 appears to rebuild the projection table per record. The probe (lib/mix/tasks/hrdps_probe.ex, deleted in 1da78a80) had already proved that wgrib2 -lon at five Canadian cities completes in under a second. Add decoder::extract_points that uses the same -lon flag for the HRDPS-only point set (~57k cells), batched at 2000 points per wgrib2 invocation to stay under Linux ARG_MAX. run_chain_step_hrdps now calls extract_points instead of extract_grid. The post-extract HashSet filter to hrdps_only_points is gone — extraction is already restricted to those cells. Per-cell wall time should drop from "never completes" to ~30-90 s for 57k cells. extract_grid stays unchanged — it works correctly for HRRR's CONUS Lambert source. Tests: parse_lon_output exercises the wgrib2 -s -lon output shape and proves multi-record / multi-point parsing. parse_lon_segment unit tests cover the undefined-value sentinel and the snap-back rejection that ignores cells wgrib2 returned that don't correspond to anything in the requested batch. --- rust/prop_grid_rs/src/decoder.rs | 222 ++++++++++++++++++++++++++++++ rust/prop_grid_rs/src/pipeline.rs | 33 ++--- 2 files changed, 237 insertions(+), 18 deletions(-) diff --git a/rust/prop_grid_rs/src/decoder.rs b/rust/prop_grid_rs/src/decoder.rs index 358c514f..f8f1cbf1 100644 --- a/rust/prop_grid_rs/src/decoder.rs +++ b/rust/prop_grid_rs/src/decoder.rs @@ -88,6 +88,190 @@ pub fn denormalize_lon(lon: f64) -> f64 { } } +/// Per-cell HRDPS extraction via `wgrib2 -lon LON LAT`. +/// +/// `extract_grid` (`-lola`) interpolates the entire source grid onto a +/// target rectangle, which on HRDPS's rotated lat/lon source runs at +/// 10+ minutes per chain step (production observation 2026-04-29). The +/// `-lon` per-point path that the probe used is far cheaper because it +/// only computes the rotation math at the points we ask for. For +/// 57k Canadian-only cells this drops chain-step wall time from "never +/// completes" to roughly 30-90 seconds. +/// +/// Batched at `POINT_BATCH` points per wgrib2 invocation to keep the +/// argv list comfortably under Linux's `ARG_MAX` (about 2 MiB) and keep +/// stdout parseable in chunks. +pub fn extract_points( + grib: &[u8], + match_pattern: &str, + points: &[(f64, f64)], +) -> Result { + let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?; + let nanos = SystemTime::now() + .duration_since(UNIX_EPOCH) + .map(|d| d.as_nanos()) + .unwrap_or(0); + let uniq = UNIQUE.fetch_add(1, Ordering::Relaxed); + let grib_path: PathBuf = std::env::temp_dir().join(format!("hrdps_{nanos}_{uniq}.grib2")); + + std::fs::write(&grib_path, grib)?; + let result = run_wgrib2_lon_batched(&wgrib2, &grib_path, match_pattern, points); + let _ = std::fs::remove_file(&grib_path); + result +} + +const POINT_BATCH: usize = 2_000; + +fn run_wgrib2_lon_batched( + wgrib2: &Path, + grib_path: &Path, + match_pattern: &str, + points: &[(f64, f64)], +) -> Result { + let mut out: PointGrid = HashMap::with_capacity(points.len()); + + for batch in points.chunks(POINT_BATCH) { + let batch_grid = run_wgrib2_lon_one(wgrib2, grib_path, match_pattern, batch)?; + for (key, cell) in batch_grid { + out.entry(key).or_default().extend(cell); + } + } + Ok(out) +} + +fn run_wgrib2_lon_one( + wgrib2: &Path, + grib_path: &Path, + match_pattern: &str, + batch: &[(f64, f64)], +) -> Result { + let mut cmd = Command::new(wgrib2); + cmd.arg(grib_path).arg("-s").arg("-match").arg(match_pattern); + for &(lat, lon) in batch { + cmd.arg("-lon").arg(format!("{lon}")).arg(format!("{lat}")); + } + + let Output { + status, + stdout, + stderr, + } = cmd.output()?; + + if !status.success() { + let mut combined = stdout; + combined.extend_from_slice(&stderr); + let text = String::from_utf8_lossy(&combined).to_string(); + let snippet: String = text.chars().take(200).collect(); + return Err(DecodeError::Wgrib2Failed { + code: status.code().unwrap_or(-1), + stderr: snippet, + }); + } + + let text = String::from_utf8_lossy(&stdout); + Ok(parse_lon_output(&text, batch)) +} + +/// Parse `wgrib2 -s -lon` text output. Each record line is a colon- +/// separated header followed by per-point `lon=…,lat=…,val=…` segments, +/// one per `-lon` arg. +/// +/// Example line: +/// `1:0:d=...:TMP:2 m above ground:anl:lon=280.369,lat=43.670,val=280.97:lon=...` +fn parse_lon_output(text: &str, batch: &[(f64, f64)]) -> PointGrid { + let mut out: PointGrid = HashMap::new(); + + for line in text.lines() { + if line.is_empty() { + continue; + } + let mut parts = line.split(':'); + let _n = parts.next(); + let _offset = parts.next(); + let _date = parts.next(); + let var = match parts.next() { + Some(v) => v.to_string(), + None => continue, + }; + let level = match parts.next() { + Some(l) => l.to_string(), + None => continue, + }; + // Skip the record-type/forecast-time field (anl, "240 min fcst", …). + let _kind = parts.next(); + let key: std::sync::Arc = format!("{var}:{level}").into(); + + for segment in parts { + if let Some(value) = parse_lon_segment(segment, batch) { + let (lat_key, lon_key, val) = value; + out.entry((lat_key, lon_key)) + .or_default() + .insert(std::sync::Arc::clone(&key), val); + } + } + } + + out +} + +/// `lon=242.958,lat=32.938,val=306.5` → (lat_key, lon_key, val). +/// Snaps to the nearest point in `batch` (handles wgrib2's lon-180/360 +/// convention by checking both signs). +fn parse_lon_segment(segment: &str, batch: &[(f64, f64)]) -> Option<(i32, i32, f32)> { + let mut lon_part = None; + let mut lat_part = None; + let mut val_part = None; + for kv in segment.split(',') { + let (k, v) = kv.split_once('=')?; + let parsed: f64 = v.trim().parse().ok()?; + match k.trim() { + "lon" => lon_part = Some(parsed), + "lat" => lat_part = Some(parsed), + "val" => val_part = Some(parsed), + _ => {} + } + } + let raw_lon = lon_part?; + let raw_lat = lat_part?; + let val = val_part?; + + if val > UNDEFINED_VALUE as f64 / 2.0 { + return None; + } + + let target_lon = denormalize_lon(raw_lon); + let snapped = snap_to_batch(raw_lat, target_lon, batch)?; + let lat_key = (snapped.0 * 1000.0).round() as i32; + let lon_key = (snapped.1 * 1000.0).round() as i32; + Some((lat_key, lon_key, val as f32)) +} + +fn snap_to_batch(lat: f64, lon: f64, batch: &[(f64, f64)]) -> Option<(f64, f64)> { + // wgrib2 -lon snaps to the nearest source-grid cell, so the returned + // lat/lon is close to but not exactly equal to what we asked for. + // Match it back to our request (which is on the 0.125° propagation + // grid) via the closest batch entry. + let mut best: Option<((f64, f64), f64)> = None; + for &(p_lat, p_lon) in batch { + let dlat = p_lat - lat; + let dlon = p_lon - lon; + let d2 = dlat * dlat + dlon * dlon; + match best { + None => best = Some(((p_lat, p_lon), d2)), + Some((_, b)) if d2 < b => best = Some(((p_lat, p_lon), d2)), + _ => {} + } + } + let (winner, d2) = best?; + // Reject anything more than 0.5° away — wgrib2 returned a cell that + // doesn't correspond to a point we asked for (probably a rounding + // edge case at the bbox boundary). + if d2 > 0.25 { + return None; + } + Some(winner) +} + /// Extract values at the given grid spec from an in-memory GRIB2 blob. /// Writes the blob to a temp file, runs wgrib2, parses the output. pub fn extract_grid( @@ -293,6 +477,44 @@ mod tests { assert_eq!(msgs[0].var, "TMP"); } + #[test] + fn parse_lon_output_pulls_var_level_and_snaps_to_batch() { + // Match the shape `wgrib2 -s -lon` produces. Two records, two + // points each. Note: lon=280.369 is wgrib2's 0-360 form for + // -79.631; parse_lon_segment denormalizes it. + let text = "1:0:d=2026042912:TMP:2 m above ground:anl:\ + lon=280.369,lat=43.670,val=280.97:\ + lon=235.000,lat=49.190,val=281.62\n\ + 2:3527800:d=2026042912:DPT:2 m above ground:anl:\ + lon=280.369,lat=43.670,val=275.00:\ + lon=235.000,lat=49.190,val=270.00\n"; + let batch = vec![(43.670, -79.631), (49.190, -125.000)]; + + let grid = parse_lon_output(text, &batch); + + let toronto_key = ((43.670_f64 * 1000.0).round() as i32, (-79.631_f64 * 1000.0).round() as i32); + let toronto = grid.get(&toronto_key).expect("toronto cell"); + let tmp = toronto.get("TMP:2 m above ground").copied().expect("tmp present"); + assert!((tmp - 280.97).abs() < 0.01); + let dpt = toronto.get("DPT:2 m above ground").copied().expect("dpt present"); + assert!((dpt - 275.00).abs() < 0.01); + } + + #[test] + fn parse_lon_segment_drops_undefined_values() { + let batch = vec![(43.670, -79.631)]; + let segment = "lon=280.369,lat=43.670,val=9.999e20"; + assert!(parse_lon_segment(segment, &batch).is_none()); + } + + #[test] + fn parse_lon_segment_rejects_far_off_points() { + let batch = vec![(43.670, -79.631)]; + // wgrib2 returned a cell for somewhere completely unrelated. + let segment = "lon=240.0,lat=10.0,val=290.0"; + assert!(parse_lon_segment(segment, &batch).is_none()); + } + #[test] fn parse_lola_binary_reconstructs_grid() { // Build a synthetic 2-message, 3×2 grid. Layout matches wgrib2: diff --git a/rust/prop_grid_rs/src/pipeline.rs b/rust/prop_grid_rs/src/pipeline.rs index 8cd1c79c..4dd0b7be 100644 --- a/rust/prop_grid_rs/src/pipeline.rs +++ b/rust/prop_grid_rs/src/pipeline.rs @@ -15,7 +15,7 @@ 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_grid_spec, hrdps_only_points, wgrib2_grid_spec}; +use crate::grid::{hrdps_only_points, wgrib2_grid_spec}; use crate::hrdps_fetcher::{self, HrdpsClient}; use crate::native_duct; use crate::nexrad::{self, NexradObservation}; @@ -305,7 +305,6 @@ pub async fn run_chain_step_hrdps( 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 @@ -313,30 +312,28 @@ pub async fn run_chain_step_hrdps( // 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_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) + let grid = tokio::task::spawn_blocking(move || { + decoder::extract_points(&blob_for_decode, &pattern_owned, &canadian_points) }) .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; + 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());