Graham McIntire
db57b3a1c7
Item 1 — .sgrid dense binary scalar format replacing chunked gzip+msgpack - New rust/prop_grid_rs/src/sgrid.rs: magic SGRD, same header layout as .pgrid, 20-field cell-major f32 body, write_atomic (tmp+rename, NFS-safe) - New lib/microwaveprop/weather/sgrid.ex: Elixir reader modeled on pgrid.ex (pread single-cell reads, bounds-filtered viewport reads, NaN→nil sentinel) - ScalarFile updated to prefer .sgrid reads, chunked .mp.gz as fallback - Pipeline writes .sgrid alongside existing chunked format (all three paths) Item 2 — HRDPS decode-once + rotated-pole index, restore 0.125° resolution - New rust/prop_grid_rs/src/rotated_pole.rs: CF-convention geographic→rotated transform, OnceLock-cached GDS params parsed from wgrib2 -grid, precomputed Vec<u32> lookup table mapping target cells to native grid indices - Native decode path in decoder.rs: wgrib2 -no_header -order we:sn -bin (raw f32 dump, ~0.32 s/message) + indexing via lookup table - HRDPS_STEP: 0.5° → 0.125° (4× finer, ~57k Canadian cells vs ~3.5k) Item 3 — k8s memory limits: 3Gi → 1.5Gi (per-task grid footprint: ~200-400 MB HashMap → ~18 MB dense planes) Item 4 — CLAUDE.md and profiles_file.ex documentation drift fixed: .pgrid primary, .mp.gz legacy, .sgrid added, write_atomic protocol doc, cleanup gaps reorganized, 'Only f00 is persisted' corrected to f00..f48
820 lines
28 KiB
Rust
820 lines
28 KiB
Rust
//! wgrib2 subprocess wrapper. 1:1 port of the hot paths in
|
||
//! `lib/microwaveprop/weather/grib2/wgrib2.ex` used by f01..f48:
|
||
//! * `extract_grid` (write GRIB2 to tmp → `wgrib2 -lola … bin` → parse)
|
||
//! * inventory parsing from stdout
|
||
//! * `parse_lola_binary` (Fortran-unformatted IEEE 754 little-endian f32)
|
||
//!
|
||
//! wgrib2 is called as a child process exactly as the Elixir version
|
||
//! does — same CLI flags, same temp-file lifecycle. Output binary is
|
||
//! stride-addressed: each message = 4-byte LE record length + N*4 bytes
|
||
//! of data + 4-byte duplicate record length.
|
||
|
||
use std::collections::HashMap;
|
||
use std::path::{Path, PathBuf};
|
||
use std::process::{Command, Output};
|
||
use std::sync::atomic::{AtomicU64, Ordering};
|
||
use std::time::{SystemTime, UNIX_EPOCH};
|
||
|
||
use crate::field_grid::FieldGrid;
|
||
use crate::grid::GridSpec;
|
||
|
||
const UNDEFINED_VALUE: f32 = 9.999e20;
|
||
static UNIQUE: AtomicU64 = AtomicU64::new(1);
|
||
|
||
#[derive(Debug, thiserror::Error)]
|
||
pub enum DecodeError {
|
||
#[error("wgrib2 not available on PATH")]
|
||
Wgrib2NotAvailable,
|
||
#[error("wgrib2 failed (exit {code}): {stderr}")]
|
||
Wgrib2Failed { code: i32, stderr: String },
|
||
#[error("io: {0}")]
|
||
Io(#[from] std::io::Error),
|
||
}
|
||
|
||
#[derive(Debug, Clone)]
|
||
pub struct Message {
|
||
pub var: String,
|
||
pub level: String,
|
||
}
|
||
|
||
// The former `CellValues` / `PointGrid` aliases (a HashMap of ~95 k
|
||
// per-cell HashMaps) are gone — see `field_grid::FieldGrid`, which
|
||
// stores the same data as dense per-message planes. Decoded grids are
|
||
// `FieldGrid` everywhere now.
|
||
|
||
pub fn wgrib2_available() -> bool {
|
||
which_wgrib2().is_some()
|
||
}
|
||
|
||
fn which_wgrib2() -> Option<PathBuf> {
|
||
if let Ok(val) = std::env::var("WGRIB2") {
|
||
return Some(PathBuf::from(val));
|
||
}
|
||
let path_var = std::env::var_os("PATH")?;
|
||
for dir in std::env::split_paths(&path_var) {
|
||
let candidate = dir.join("wgrib2");
|
||
if candidate.is_file() {
|
||
return Some(candidate);
|
||
}
|
||
}
|
||
None
|
||
}
|
||
|
||
pub fn normalize_lon(lon: f64) -> f64 {
|
||
if lon < 0.0 {
|
||
lon + 360.0
|
||
} else {
|
||
lon
|
||
}
|
||
}
|
||
|
||
pub fn denormalize_lon(lon: f64) -> f64 {
|
||
if lon > 180.0 {
|
||
lon - 360.0
|
||
} else {
|
||
lon
|
||
}
|
||
}
|
||
|
||
/// 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)],
|
||
grid_spec: GridSpec,
|
||
) -> Result<FieldGrid, 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, grid_spec);
|
||
let _ = std::fs::remove_file(&grib_path);
|
||
result
|
||
}
|
||
|
||
// 1000 points per wgrib2 invocation. With ~57k Canadian cells that's
|
||
// 57 invocations per chain step — enough to amortize process-spawn
|
||
// overhead, small enough that any future wgrib2 regression around
|
||
// large -lon arg lists trips on a manageable failure mode.
|
||
//
|
||
// Note: wgrib2 3.6.0 had a memory-corruption bug surfaced by HRDPS
|
||
// rotated lat/lon files at *any* batch size (verified at N=200 in
|
||
// production on 2026-04-30 — `free(): invalid size`). The fix was
|
||
// bumping the Dockerfile to wgrib2 3.8.0; this batch size is just
|
||
// a normal performance tuning choice, not a workaround.
|
||
const POINT_BATCH: usize = 1_000;
|
||
|
||
fn run_wgrib2_lon_batched(
|
||
wgrib2: &Path,
|
||
grib_path: &Path,
|
||
match_pattern: &str,
|
||
points: &[(f64, f64)],
|
||
grid_spec: GridSpec,
|
||
) -> Result<FieldGrid, DecodeError> {
|
||
let mut out = FieldGrid::new(grid_spec);
|
||
|
||
for batch in points.chunks(POINT_BATCH) {
|
||
let batch_values = run_wgrib2_lon_one(wgrib2, grib_path, match_pattern, batch, &grid_spec)?;
|
||
for (name, cells) in batch_values {
|
||
let (_, plane) = out.plane_mut_or_insert(&name);
|
||
for (cell, value) in cells {
|
||
plane[cell] = value;
|
||
}
|
||
}
|
||
}
|
||
Ok(out)
|
||
}
|
||
|
||
/// One wgrib2 `-lon` invocation, keyed `"VAR:LEVEL" -> [(cell_index, value)]`.
|
||
/// Values are gathered per-name so the caller can scatter each name into a
|
||
/// single plane instead of touching a per-cell map.
|
||
type LonBatchValues = HashMap<String, Vec<(usize, f32)>>;
|
||
|
||
fn run_wgrib2_lon_one(
|
||
wgrib2: &Path,
|
||
grib_path: &Path,
|
||
match_pattern: &str,
|
||
batch: &[(f64, f64)],
|
||
grid_spec: &GridSpec,
|
||
) -> Result<LonBatchValues, 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() {
|
||
// Surface the FIRST chars of stderr (the actual diagnostic)
|
||
// rather than stdout (which is mostly per-cell numeric output
|
||
// that drowned the interesting failure context). status.code()
|
||
// returns None on signal-kill — flag that with -1 so the
|
||
// failure mode is distinguishable from a real exit code.
|
||
let stderr_text = String::from_utf8_lossy(&stderr).to_string();
|
||
let stderr_snippet: String = stderr_text.chars().take(400).collect();
|
||
let signaled = status.code().is_none();
|
||
let suffix = if signaled {
|
||
format!(" (killed by signal, batch={})", batch.len())
|
||
} else {
|
||
String::new()
|
||
};
|
||
return Err(DecodeError::Wgrib2Failed {
|
||
code: status.code().unwrap_or(-1),
|
||
stderr: format!("{stderr_snippet}{suffix}"),
|
||
});
|
||
}
|
||
|
||
let text = String::from_utf8_lossy(&stdout);
|
||
Ok(parse_lon_output(&text, batch, grid_spec))
|
||
}
|
||
|
||
/// 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)], grid_spec: &GridSpec) -> LonBatchValues {
|
||
let mut out: LonBatchValues = 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,
|
||
None => continue,
|
||
};
|
||
let level = match parts.next() {
|
||
Some(l) => l,
|
||
None => continue,
|
||
};
|
||
// Skip the record-type/forecast-time field (anl, "240 min fcst", …).
|
||
let _kind = parts.next();
|
||
let key = format!("{var}:{level}");
|
||
let slot = out.entry(key).or_default();
|
||
|
||
for segment in parts {
|
||
if let Some((lat, lon, val)) = parse_lon_segment(segment, batch) {
|
||
if let Some(cell) = crate::field_grid::cell_index_for(grid_spec, lat, lon) {
|
||
slot.push((cell, val));
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
out
|
||
}
|
||
|
||
/// `lon=242.958,lat=32.938,val=306.5` → (lat, lon, val) snapped to the
|
||
/// nearest point in `batch` (handles wgrib2's lon-180/360 convention by
|
||
/// denormalising before matching).
|
||
fn parse_lon_segment(segment: &str, batch: &[(f64, f64)]) -> Option<(f64, f64, 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 (lat, lon) = snap_to_batch(raw_lat, target_lon, batch)?;
|
||
Some((lat, lon, 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(
|
||
grib: &[u8],
|
||
match_pattern: &str,
|
||
grid_spec: GridSpec,
|
||
) -> Result<FieldGrid, DecodeError> {
|
||
let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?;
|
||
let tmp_dir = std::env::temp_dir();
|
||
let uniq = UNIQUE.fetch_add(1, Ordering::Relaxed);
|
||
let nanos = SystemTime::now()
|
||
.duration_since(UNIX_EPOCH)
|
||
.map(|d| d.as_nanos())
|
||
.unwrap_or(0);
|
||
let grib_path = tmp_dir.join(format!("hrrr_{nanos}_{uniq}.grib2"));
|
||
let bin_path = tmp_dir.join(format!("hrrr_{nanos}_{uniq}.grib2.lola.bin"));
|
||
|
||
let result: Result<FieldGrid, DecodeError> = (|| {
|
||
std::fs::write(&grib_path, grib)?;
|
||
run_wgrib2_lola(&wgrib2, &grib_path, match_pattern, &grid_spec, &bin_path)
|
||
})();
|
||
|
||
let _ = std::fs::remove_file(&grib_path);
|
||
let _ = std::fs::remove_file(&bin_path);
|
||
result
|
||
}
|
||
|
||
/// File-path variant — avoids the intermediate `write(tmp, binary)` copy.
|
||
pub fn extract_grid_from_file(
|
||
grib_path: &Path,
|
||
match_pattern: &str,
|
||
grid_spec: GridSpec,
|
||
) -> Result<FieldGrid, DecodeError> {
|
||
let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?;
|
||
let bin_path = sibling_tmp_path(grib_path, "lola");
|
||
let result = run_wgrib2_lola(&wgrib2, grib_path, match_pattern, &grid_spec, &bin_path);
|
||
let _ = std::fs::remove_file(&bin_path);
|
||
result
|
||
}
|
||
|
||
fn sibling_tmp_path(path: &Path, suffix: &str) -> PathBuf {
|
||
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 mut s = path.as_os_str().to_owned();
|
||
s.push(format!(".{suffix}.{nanos}.{uniq}.bin"));
|
||
PathBuf::from(s)
|
||
}
|
||
|
||
fn run_wgrib2_lola(
|
||
wgrib2: &Path,
|
||
grib_path: &Path,
|
||
match_pattern: &str,
|
||
grid_spec: &GridSpec,
|
||
bin_path: &Path,
|
||
) -> Result<FieldGrid, DecodeError> {
|
||
let lon_spec = format!(
|
||
"{}:{}:{}",
|
||
normalize_lon(grid_spec.lon_start),
|
||
grid_spec.lon_count,
|
||
grid_spec.lon_step
|
||
);
|
||
let lat_spec = format!(
|
||
"{}:{}:{}",
|
||
grid_spec.lat_start, grid_spec.lat_count, grid_spec.lat_step
|
||
);
|
||
|
||
let Output {
|
||
status,
|
||
stdout,
|
||
stderr,
|
||
} = Command::new(wgrib2)
|
||
.arg(grib_path)
|
||
.arg("-match")
|
||
.arg(match_pattern)
|
||
.arg("-lola")
|
||
.arg(&lon_spec)
|
||
.arg(&lat_spec)
|
||
.arg(bin_path)
|
||
.arg("bin")
|
||
.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 inventory = String::from_utf8_lossy(&stdout);
|
||
let messages = parse_inventory(&inventory);
|
||
|
||
match std::fs::read(bin_path) {
|
||
Ok(bin) => Ok(parse_lola_binary(&bin, &messages, grid_spec)),
|
||
// wgrib2 matched nothing, so it never created the output file.
|
||
// An empty grid (right spec, zero planes) is the correct result.
|
||
Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(FieldGrid::new(*grid_spec)),
|
||
Err(e) => Err(e.into()),
|
||
}
|
||
}
|
||
|
||
/// Parse `wgrib2`'s stdout inventory (`"n:offset:date:var:level:..."`).
|
||
/// Malformed lines are dropped.
|
||
pub fn parse_inventory(out: &str) -> Vec<Message> {
|
||
out.lines()
|
||
.filter_map(|line| {
|
||
let mut parts = line.splitn(8, ':');
|
||
let _n = parts.next()?;
|
||
let _offset = parts.next()?;
|
||
let _date = parts.next()?;
|
||
let var = parts.next()?;
|
||
let level = parts.next()?;
|
||
Some(Message {
|
||
var: var.to_string(),
|
||
level: level.to_string(),
|
||
})
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
/// Parse the Fortran-unformatted output file from `wgrib2 -lola`.
|
||
/// Layout: per message = 4-byte LE u32 record length + N·4 bytes f32 LE
|
||
/// data + 4-byte LE u32 record length. `N = nx * ny`. Cells whose value
|
||
/// exceeds `UNDEFINED_VALUE / 2` are replaced with `NaN` (the sentinel
|
||
/// wgrib2 writes for "no data"; `FieldGrid::at` maps NaN back to `None`,
|
||
/// matching the old behaviour of omitting the key entirely).
|
||
///
|
||
/// wgrib2 writes each record row-major over (lat, lon) in exactly the
|
||
/// order `FieldGrid` indexes cells, so the body is a bulk
|
||
/// little-endian `f32` decode per message — no per-cell coordinate
|
||
/// arithmetic, no hashing, no allocation beyond the plane itself.
|
||
pub fn parse_lola_binary(bin: &[u8], messages: &[Message], grid_spec: &GridSpec) -> FieldGrid {
|
||
let nx = grid_spec.lon_count;
|
||
let ny = grid_spec.lat_count;
|
||
let n_cells = nx * ny;
|
||
let bytes_per_msg = n_cells * 4;
|
||
let record_overhead = 8;
|
||
let stride = bytes_per_msg + record_overhead;
|
||
|
||
let mut out = FieldGrid::new(*grid_spec);
|
||
|
||
for (msg_idx, msg) in messages.iter().enumerate() {
|
||
let data_offset = msg_idx * stride + 4;
|
||
if data_offset + bytes_per_msg > bin.len() {
|
||
continue;
|
||
}
|
||
let chunk = &bin[data_offset..data_offset + bytes_per_msg];
|
||
let mut plane = Vec::with_capacity(n_cells);
|
||
plane.extend(chunk.chunks_exact(4).map(|b| {
|
||
// try_into on a chunks_exact(4) window is statically length 4.
|
||
let v = f32::from_le_bytes(b.try_into().unwrap());
|
||
if v > UNDEFINED_VALUE / 2.0 {
|
||
f32::NAN
|
||
} else {
|
||
v
|
||
}
|
||
}));
|
||
out.push_plane(&format!("{}:{}", msg.var, msg.level), plane);
|
||
}
|
||
|
||
out
|
||
}
|
||
|
||
/// Decode an HRDPS rotated-pole GRIB2 blob using decode-once +
|
||
/// lookup-many. Replaces the per-point `-lon` extraction that used to
|
||
/// brute-force the rotation math once per output cell.
|
||
///
|
||
/// 1. Runs `wgrib2 -no_header -order we:sn -bin` to dump the native
|
||
/// grid as raw f32 blocks (~0.32 s per message).
|
||
/// 2. Parses the inventory from stdout to get variable+level names.
|
||
/// 3. For each message, builds a dense `FieldGrid` plane by indexing
|
||
/// into the native grid via the precomputed `lookup` table.
|
||
///
|
||
/// `native_nx` × `native_ny` is the native grid size (2540×1290 for
|
||
/// HRDPS 0.0225°). `lookup` maps each cell in `target_spec` to a
|
||
/// native grid index, or `None` if the target cell is outside the
|
||
/// native domain.
|
||
pub fn decode_hrdps_native(
|
||
grib: &[u8],
|
||
match_pattern: &str,
|
||
target_spec: &GridSpec,
|
||
lookup: &[Option<u32>],
|
||
native_nx: u32,
|
||
native_ny: u32,
|
||
) -> Result<FieldGrid, DecodeError> {
|
||
let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?;
|
||
let tmp_dir = std::env::temp_dir();
|
||
let uniq = UNIQUE.fetch_add(1, Ordering::Relaxed);
|
||
let nanos = SystemTime::now()
|
||
.duration_since(UNIX_EPOCH)
|
||
.map(|d| d.as_nanos())
|
||
.unwrap_or(0);
|
||
let grib_path = tmp_dir.join(format!("hrdps_native_{nanos}_{uniq}.grib2"));
|
||
let bin_path = tmp_dir.join(format!("hrdps_native_{nanos}_{uniq}.bin"));
|
||
|
||
std::fs::write(&grib_path, grib)?;
|
||
|
||
let result = run_wgrib2_native_bin(
|
||
&wgrib2,
|
||
&grib_path,
|
||
match_pattern,
|
||
&bin_path,
|
||
target_spec,
|
||
lookup,
|
||
(native_nx, native_ny),
|
||
);
|
||
|
||
let _ = std::fs::remove_file(&grib_path);
|
||
let _ = std::fs::remove_file(&bin_path);
|
||
result
|
||
}
|
||
|
||
fn run_wgrib2_native_bin(
|
||
wgrib2: &Path,
|
||
grib_path: &Path,
|
||
match_pattern: &str,
|
||
bin_path: &Path,
|
||
target_spec: &GridSpec,
|
||
lookup: &[Option<u32>],
|
||
native_dims: (u32, u32),
|
||
) -> Result<FieldGrid, DecodeError> {
|
||
let Output {
|
||
status,
|
||
stdout,
|
||
stderr,
|
||
} = Command::new(wgrib2)
|
||
.arg(grib_path)
|
||
.arg("-match")
|
||
.arg(match_pattern)
|
||
.arg("-no_header")
|
||
.arg("-order")
|
||
.arg("we:sn")
|
||
.arg("-bin")
|
||
.arg(bin_path)
|
||
.output()?;
|
||
|
||
if !status.success() {
|
||
let stderr_text = String::from_utf8_lossy(&stderr).to_string();
|
||
let snippet: String = stderr_text.chars().take(400).collect();
|
||
return Err(DecodeError::Wgrib2Failed {
|
||
code: status.code().unwrap_or(-1),
|
||
stderr: snippet,
|
||
});
|
||
}
|
||
|
||
let inventory = String::from_utf8_lossy(&stdout);
|
||
let messages = parse_inventory(&inventory);
|
||
|
||
match std::fs::read(bin_path) {
|
||
Ok(bin) => Ok(parse_native_binary(
|
||
&bin,
|
||
&messages,
|
||
target_spec,
|
||
lookup,
|
||
native_dims,
|
||
)),
|
||
Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(FieldGrid::new(*target_spec)),
|
||
Err(e) => Err(e.into()),
|
||
}
|
||
}
|
||
|
||
/// Parse the raw `-no_header -order we:sn -bin` output: each message is
|
||
/// `nx * ny * 4` bytes of f32 LE, one block per message in the same
|
||
/// order as the inventory. No record markers, no headers.
|
||
fn parse_native_binary(
|
||
bin: &[u8],
|
||
messages: &[Message],
|
||
target_spec: &GridSpec,
|
||
lookup: &[Option<u32>],
|
||
(native_nx, native_ny): (u32, u32),
|
||
) -> FieldGrid {
|
||
let n_cells = target_spec.lat_count * target_spec.lon_count;
|
||
let bytes_per_msg = (native_nx as usize) * (native_ny as usize) * 4;
|
||
let n_native = (native_nx as usize) * (native_ny as usize);
|
||
|
||
let mut out = FieldGrid::new(*target_spec);
|
||
|
||
for (msg_idx, msg) in messages.iter().enumerate() {
|
||
let data_offset = msg_idx * bytes_per_msg;
|
||
if data_offset + bytes_per_msg > bin.len() {
|
||
continue;
|
||
}
|
||
let chunk = &bin[data_offset..data_offset + bytes_per_msg];
|
||
|
||
// Decode the entire native grid into f32 values once per message.
|
||
// For the 2540×1290 (= 3.28M cell) native grid, this is ~13 MB
|
||
// per message, which is acceptable for the ~41 messages HRDPS
|
||
// publishes. The alternative would be indexing into raw bytes per
|
||
// target cell, which adds per-cell bounds checks and offset math.
|
||
let native: Vec<f32> = chunk
|
||
.chunks_exact(4)
|
||
.map(|b| {
|
||
let v = f32::from_le_bytes(b.try_into().unwrap());
|
||
if v > UNDEFINED_VALUE / 2.0 {
|
||
f32::NAN
|
||
} else {
|
||
v
|
||
}
|
||
})
|
||
.collect();
|
||
debug_assert_eq!(native.len(), n_native);
|
||
|
||
// Build the target plane by indexing into the native grid.
|
||
let name = format!("{}:{}", msg.var, msg.level);
|
||
let mut plane = vec![f32::NAN; n_cells];
|
||
for (cell, slot) in lookup.iter().enumerate() {
|
||
if let Some(native_cell) = slot {
|
||
if (*native_cell as usize) < n_native {
|
||
plane[cell] = native[*native_cell as usize];
|
||
}
|
||
}
|
||
}
|
||
out.push_plane(&name, plane);
|
||
}
|
||
|
||
out
|
||
}
|
||
|
||
pub fn key_to_latlon(key: (i32, i32)) -> (f64, f64) {
|
||
(key.0 as f64 / 1000.0, key.1 as f64 / 1000.0)
|
||
}
|
||
|
||
/// Run `wgrib2 -grid` on a GRIB2 file and return the stdout text.
|
||
/// Used to extract rotated-pole grid definition parameters at startup
|
||
/// so the lookup table can be validated against the actual GDS.
|
||
pub fn extract_grid_metadata(grib_path: &Path) -> Result<String, DecodeError> {
|
||
let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?;
|
||
let Output {
|
||
status,
|
||
stdout,
|
||
stderr,
|
||
} = Command::new(&wgrib2).arg(grib_path).arg("-grid").output()?;
|
||
|
||
if !status.success() {
|
||
let stderr_text = String::from_utf8_lossy(&stderr).to_string();
|
||
let snippet: String = stderr_text.chars().take(400).collect();
|
||
return Err(DecodeError::Wgrib2Failed {
|
||
code: status.code().unwrap_or(-1),
|
||
stderr: snippet,
|
||
});
|
||
}
|
||
|
||
Ok(String::from_utf8_lossy(&stdout).to_string())
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
#[test]
|
||
fn lon_round_trip() {
|
||
for lon in [-125.0, -97.0, -66.0, 0.0, 179.5] {
|
||
assert!((denormalize_lon(normalize_lon(lon)) - lon).abs() < 1e-9);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn inventory_parses_var_level() {
|
||
let text = "1:0:d=2026041915:TMP:2 m above ground:anl:\n\
|
||
2:1234:d=2026041915:DPT:2 m above ground:anl:\n";
|
||
let msgs = parse_inventory(text);
|
||
assert_eq!(msgs.len(), 2);
|
||
assert_eq!(msgs[0].var, "TMP");
|
||
assert_eq!(msgs[1].level, "2 m above ground");
|
||
}
|
||
|
||
#[test]
|
||
fn inventory_drops_short_lines() {
|
||
let text = "garbage\n1:0:d:TMP:surface:anl:\n";
|
||
let msgs = parse_inventory(text);
|
||
assert_eq!(msgs.len(), 1);
|
||
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)];
|
||
// Spec chosen so both batch points land exactly on a cell.
|
||
let spec = GridSpec {
|
||
lon_start: -125.0,
|
||
lon_count: 400,
|
||
lon_step: 0.125,
|
||
lat_start: 43.0,
|
||
lat_count: 60,
|
||
lat_step: 0.125,
|
||
};
|
||
|
||
let values = parse_lon_output(text, &batch, &spec);
|
||
|
||
let toronto = crate::field_grid::cell_index_for(&spec, 43.670, -79.631)
|
||
.expect("toronto lands on the grid");
|
||
let tmp = values["TMP:2 m above ground"]
|
||
.iter()
|
||
.find(|(c, _)| *c == toronto)
|
||
.expect("tmp present")
|
||
.1;
|
||
assert!((tmp - 280.97).abs() < 0.01);
|
||
let dpt = values["DPT:2 m above ground"]
|
||
.iter()
|
||
.find(|(c, _)| *c == toronto)
|
||
.expect("dpt present")
|
||
.1;
|
||
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:
|
||
// [4 bytes len][6 f32 values][4 bytes len] per message.
|
||
let spec = GridSpec {
|
||
lon_start: -100.0,
|
||
lon_count: 3,
|
||
lon_step: 0.5,
|
||
lat_start: 30.0,
|
||
lat_count: 2,
|
||
lat_step: 0.5,
|
||
};
|
||
let msgs = vec![
|
||
Message {
|
||
var: "TMP".into(),
|
||
level: "surface".into(),
|
||
},
|
||
Message {
|
||
var: "DPT".into(),
|
||
level: "surface".into(),
|
||
},
|
||
];
|
||
let mut buf = Vec::new();
|
||
for (msg_idx, _msg) in msgs.iter().enumerate() {
|
||
let len: u32 = 3 * 2 * 4;
|
||
buf.extend_from_slice(&len.to_le_bytes());
|
||
for j in 0..2 {
|
||
for i in 0..3 {
|
||
let v = (msg_idx as f32) * 100.0 + j as f32 * 10.0 + i as f32;
|
||
buf.extend_from_slice(&v.to_le_bytes());
|
||
}
|
||
}
|
||
buf.extend_from_slice(&len.to_le_bytes());
|
||
}
|
||
|
||
let grid = parse_lola_binary(&buf, &msgs, &spec);
|
||
assert_eq!(grid.n_cells(), 6);
|
||
assert_eq!(grid.n_planes(), 2);
|
||
let tmp = grid.plane_id("TMP:surface").unwrap();
|
||
let dpt = grid.plane_id("DPT:surface").unwrap();
|
||
|
||
// Cell (30.0, -100.0) → i=0, j=0, msg0 value = 0.0, msg1 = 100.0
|
||
let k = grid.cell_index(30.0, -100.0).unwrap();
|
||
assert_eq!(grid.at(tmp, k), Some(0.0));
|
||
assert_eq!(grid.at(dpt, k), Some(100.0));
|
||
|
||
// Cell (30.5, -99.0) → i=2, j=1, msg0 = 12.0
|
||
let k2 = grid.cell_index(30.5, -99.0).unwrap();
|
||
assert_eq!(grid.at(tmp, k2), Some(12.0));
|
||
assert_eq!(grid.at(dpt, k2), Some(112.0));
|
||
}
|
||
|
||
#[test]
|
||
fn parse_lola_skips_undefined_sentinel() {
|
||
let spec = GridSpec {
|
||
lon_start: -100.0,
|
||
lon_count: 2,
|
||
lon_step: 1.0,
|
||
lat_start: 30.0,
|
||
lat_count: 1,
|
||
lat_step: 1.0,
|
||
};
|
||
let msgs = vec![Message {
|
||
var: "TMP".into(),
|
||
level: "surface".into(),
|
||
}];
|
||
let mut buf = Vec::new();
|
||
let len: u32 = 2 * 4;
|
||
buf.extend_from_slice(&len.to_le_bytes());
|
||
buf.extend_from_slice(&UNDEFINED_VALUE.to_le_bytes()); // i=0
|
||
buf.extend_from_slice(&42.0f32.to_le_bytes()); // i=1
|
||
buf.extend_from_slice(&len.to_le_bytes());
|
||
|
||
let grid = parse_lola_binary(&buf, &msgs, &spec);
|
||
let tmp = grid.plane_id("TMP:surface").unwrap();
|
||
// i=0 carried the sentinel and must read back as missing; only
|
||
// i=1 survived. The old representation expressed this by never
|
||
// inserting the key for i=0.
|
||
assert_eq!(grid.at(tmp, grid.cell_index(30.0, -100.0).unwrap()), None);
|
||
assert_eq!(
|
||
grid.at(tmp, grid.cell_index(30.0, -99.0).unwrap()),
|
||
Some(42.0)
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn key_to_latlon_roundtrip() {
|
||
let (lat, lon) = key_to_latlon((32_500, -97_250));
|
||
assert_eq!(lat, 32.5);
|
||
assert_eq!(lon, -97.25);
|
||
}
|
||
}
|