fix(rust/hrdps): replace -lola full-grid extract with -lon point extract
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.
This commit is contained in:
parent
5e9513bf02
commit
2348c48c26
2 changed files with 237 additions and 18 deletions
|
|
@ -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<PointGrid, DecodeError> {
|
||||
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<PointGrid, DecodeError> {
|
||||
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<PointGrid, DecodeError> {
|
||||
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<str> = 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:
|
||||
|
|
|
|||
|
|
@ -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());
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue