fix(hrdps): parse wgrib2 grid_template=1 rotated-pole output
Every HRDPS chain step failed at decode with "could not parse nx×ny
from wgrib2 -grid output", so no .hrdps.prop or .hrdps.sgrid artifact
was ever produced and the Canadian half of the weather overlay was
empty.
parse_wgrib2_grid only handled the verbose grid_template=10 shape
(Rotated Lat/lon Grid:, LatFirst:, Di:, South Pole Location:). The
wgrib2 3.8.0 binary in the pipeline image actually emits the compact
grid_template=1 shape for MSC HRDPS files:
rotated lat-lon grid:(2540 x 1290) units 1e-06 ...
lat -12.302501 to 16.700001 by 0.022500
lon 345.178780 to 42.306283 by 0.022500 #points=3276600
south pole lat=-36.088520 lon=245.305142 angle of rot=0.000000
Accept both, keeping the strict-by-default behaviour so a real MSC grid
change still surfaces as an error. Anisotropic grids are rejected
outright because RotatedPoleParams carries a single isotropic step that
geo_to_rotated applies to both axes.
Test uses the verbatim -grid output captured from the production pod.
This commit is contained in:
parent
e1585cbdcb
commit
3d214f60f0
1 changed files with 167 additions and 38 deletions
|
|
@ -1,8 +1,9 @@
|
||||||
//! Rotated lat/lon → geographic coordinate transform for HRDPS.
|
//! Rotated lat/lon → geographic coordinate transform for HRDPS.
|
||||||
//!
|
//!
|
||||||
//! HRDPS uses a rotated-pole grid (GRIB2 template 10) where the
|
//! HRDPS uses a rotated-pole grid (GRIB2 grid template 1; wgrib2 also
|
||||||
//! computational grid is regular lat/lon but the geographic grid is
|
//! renders the equivalent template 10) where the computational grid is
|
||||||
//! rotated. `wgrib2 -lola` handles this by brute-forcing per-output-point,
|
//! regular lat/lon but the geographic grid is rotated.
|
||||||
|
//! `wgrib2 -lola` handles this by brute-forcing per-output-point,
|
||||||
//! which is correct but O(output_points)-slow for the entire Canadian
|
//! which is correct but O(output_points)-slow for the entire Canadian
|
||||||
//! bbox at 0.125°.
|
//! bbox at 0.125°.
|
||||||
//!
|
//!
|
||||||
|
|
@ -175,7 +176,19 @@ fn build_lookup(params: &RotatedPoleParams, target_spec: &GridSpec) -> Vec<Optio
|
||||||
|
|
||||||
/// Parse `RotatedPoleParams` from the text output of `wgrib2 -grid`.
|
/// Parse `RotatedPoleParams` from the text output of `wgrib2 -grid`.
|
||||||
///
|
///
|
||||||
/// The output for a rotated lat/lon grid looks like (3.8.0):
|
/// wgrib2 3.8.0 prints one of two shapes for a rotated lat/lon grid
|
||||||
|
/// depending on which GDS template the file carries. Both are accepted.
|
||||||
|
///
|
||||||
|
/// `grid_template=1` — what MSC's HRDPS files actually carry:
|
||||||
|
/// ```text
|
||||||
|
/// 1:0:grid_template=1:winds(grid):
|
||||||
|
/// rotated lat-lon grid:(2540 x 1290) units 1e-06 input WE:SN output WE:SN res 56
|
||||||
|
/// lat -12.302501 to 16.700001 by 0.022500
|
||||||
|
/// lon 345.178780 to 42.306283 by 0.022500 #points=3276600
|
||||||
|
/// south pole lat=-36.088520 lon=245.305142 angle of rot=0.000000
|
||||||
|
/// ```
|
||||||
|
///
|
||||||
|
/// `grid_template=10`:
|
||||||
/// ```text
|
/// ```text
|
||||||
/// 1:0:grid_template=10:winds(N/S):
|
/// 1:0:grid_template=10:winds(N/S):
|
||||||
/// Rotated Lat/lon Grid: (2540 x 1290)
|
/// Rotated Lat/lon Grid: (2540 x 1290)
|
||||||
|
|
@ -191,46 +204,32 @@ fn build_lookup(params: &RotatedPoleParams, target_spec: &GridSpec) -> Vec<Optio
|
||||||
/// so an MSC grid change surfaces as a failure rather than silently
|
/// so an MSC grid change surfaces as a failure rather than silently
|
||||||
/// mis-indexing.
|
/// mis-indexing.
|
||||||
fn parse_wgrib2_grid(output: &str) -> Result<RotatedPoleParams, String> {
|
fn parse_wgrib2_grid(output: &str) -> Result<RotatedPoleParams, String> {
|
||||||
// Extract nx, ny from "Rotated Lat/lon Grid: (NNNN x NNNN)"
|
parse_template_10(output)
|
||||||
let (nx, ny) = output
|
.or_else(|| parse_template_1(output))
|
||||||
.lines()
|
.ok_or_else(|| {
|
||||||
.find_map(|line| {
|
format!("could not parse rotated-pole grid from wgrib2 -grid output:\n{output}")
|
||||||
let line = line.trim();
|
|
||||||
if line.starts_with("Rotated Lat/lon Grid:") {
|
|
||||||
let rest = line.strip_prefix("Rotated Lat/lon Grid:")?.trim();
|
|
||||||
let rest = rest.strip_prefix('(')?.strip_suffix(')')?;
|
|
||||||
let mut parts = rest.split('x');
|
|
||||||
let nx: u32 = parts.next()?.trim().parse().ok()?;
|
|
||||||
let ny: u32 = parts.next()?.trim().parse().ok()?;
|
|
||||||
Some((nx, ny))
|
|
||||||
} else {
|
|
||||||
None
|
|
||||||
}
|
|
||||||
})
|
})
|
||||||
.ok_or_else(|| format!("could not parse nx×ny from wgrib2 -grid output:\n{output}"))?;
|
}
|
||||||
|
|
||||||
// Extract lat0 (LatFirst), lon0 (LonFirst), d (Di)
|
/// Verbose `grid_template=10` form: one labelled `key: value` per line.
|
||||||
let lat0 = parse_float_after(output, "LatFirst:").ok_or("LatFirst not found")?;
|
fn parse_template_10(output: &str) -> Option<RotatedPoleParams> {
|
||||||
let lon0 = parse_float_after(output, "LonFirst:").ok_or("LonFirst not found")?;
|
let (nx, ny) = output.lines().find_map(|line| {
|
||||||
let d = parse_float_after(output, "Di:").ok_or("Di not found")?;
|
let rest = line.trim().strip_prefix("Rotated Lat/lon Grid:")?;
|
||||||
|
parse_dims(rest.trim())
|
||||||
|
})?;
|
||||||
|
|
||||||
// Extract south pole location
|
let lat0 = parse_float_after(output, "LatFirst:")?;
|
||||||
let sp_lat = output
|
let lon0 = parse_float_after(output, "LonFirst:")?;
|
||||||
|
let d = parse_float_after(output, "Di:")?;
|
||||||
|
|
||||||
|
let pole_line = output
|
||||||
.lines()
|
.lines()
|
||||||
.skip_while(|l| !l.trim().starts_with("South Pole Location:"))
|
.skip_while(|l| !l.trim().starts_with("South Pole Location:"))
|
||||||
.nth(1)
|
.nth(1)?;
|
||||||
.and_then(|l| parse_float_after(l, "lat:"))
|
|
||||||
.ok_or("South Pole lat not found")?;
|
|
||||||
let sp_lon = output
|
|
||||||
.lines()
|
|
||||||
.skip_while(|l| !l.trim().starts_with("South Pole Location:"))
|
|
||||||
.nth(1)
|
|
||||||
.and_then(|l| parse_float_after(l, "lon:"))
|
|
||||||
.ok_or("South Pole lon not found")?;
|
|
||||||
|
|
||||||
Ok(RotatedPoleParams {
|
Some(RotatedPoleParams {
|
||||||
sp_lat,
|
sp_lat: parse_float_after(pole_line, "lat:")?,
|
||||||
sp_lon,
|
sp_lon: parse_float_after(pole_line, "lon:")?,
|
||||||
lat0,
|
lat0,
|
||||||
lon0,
|
lon0,
|
||||||
d,
|
d,
|
||||||
|
|
@ -239,6 +238,64 @@ fn parse_wgrib2_grid(output: &str) -> Result<RotatedPoleParams, String> {
|
||||||
})
|
})
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/// Compact `grid_template=1` form. Each axis line carries the origin and
|
||||||
|
/// the step together (`lat <first> to <last> by <step>`), and the pole
|
||||||
|
/// sits on a single `key=value` line.
|
||||||
|
fn parse_template_1(output: &str) -> Option<RotatedPoleParams> {
|
||||||
|
let (nx, ny) = output.lines().find_map(|line| {
|
||||||
|
let rest = line.trim().strip_prefix("rotated lat-lon grid:")?;
|
||||||
|
parse_dims(rest.trim())
|
||||||
|
})?;
|
||||||
|
|
||||||
|
let (lat0, dj) = output.lines().find_map(|l| parse_axis(l.trim(), "lat "))?;
|
||||||
|
let (lon0, di) = output.lines().find_map(|l| parse_axis(l.trim(), "lon "))?;
|
||||||
|
|
||||||
|
// `RotatedPoleParams` carries one isotropic step because HRDPS is
|
||||||
|
// 0.0225° in both directions, and `geo_to_rotated` divides both the
|
||||||
|
// row and column offset by it. Refuse an anisotropic grid rather
|
||||||
|
// than index rows with the column step.
|
||||||
|
if (di - dj).abs() > 1e-9 {
|
||||||
|
return None;
|
||||||
|
}
|
||||||
|
|
||||||
|
let (sp_lat, sp_lon) = output.lines().find_map(|l| {
|
||||||
|
let rest = l.trim().strip_prefix("south pole ")?;
|
||||||
|
Some((
|
||||||
|
parse_float_after(rest, "lat=")?,
|
||||||
|
parse_float_after(rest, "lon=")?,
|
||||||
|
))
|
||||||
|
})?;
|
||||||
|
|
||||||
|
Some(RotatedPoleParams {
|
||||||
|
sp_lat,
|
||||||
|
sp_lon,
|
||||||
|
lat0,
|
||||||
|
lon0,
|
||||||
|
d: di,
|
||||||
|
nx,
|
||||||
|
ny,
|
||||||
|
})
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Extract `(nx, ny)` from a leading `(NNNN x NNNN)` group, ignoring any
|
||||||
|
/// trailing text on the same line.
|
||||||
|
fn parse_dims(s: &str) -> Option<(u32, u32)> {
|
||||||
|
let inner = s.strip_prefix('(')?;
|
||||||
|
let end = inner.find(')')?;
|
||||||
|
let mut parts = inner[..end].split('x');
|
||||||
|
let nx: u32 = parts.next()?.trim().parse().ok()?;
|
||||||
|
let ny: u32 = parts.next()?.trim().parse().ok()?;
|
||||||
|
Some((nx, ny))
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Parse a `<axis> <first> to <last> by <step>` line into `(first, step)`.
|
||||||
|
fn parse_axis(line: &str, prefix: &str) -> Option<(f64, f64)> {
|
||||||
|
let rest = line.strip_prefix(prefix)?;
|
||||||
|
let first: f64 = rest.split_whitespace().next()?.parse().ok()?;
|
||||||
|
let step = parse_float_after(rest, " by ")?;
|
||||||
|
Some((first, step))
|
||||||
|
}
|
||||||
|
|
||||||
fn parse_float_after(haystack: &str, prefix: &str) -> Option<f64> {
|
fn parse_float_after(haystack: &str, prefix: &str) -> Option<f64> {
|
||||||
for line in haystack.lines() {
|
for line in haystack.lines() {
|
||||||
if let Some(pos) = line.find(prefix) {
|
if let Some(pos) = line.find(prefix) {
|
||||||
|
|
@ -290,6 +347,78 @@ mod tests {
|
||||||
assert_eq!(params.ny, expected.ny);
|
assert_eq!(params.ny, expected.ny);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/// Verbatim `wgrib2 -grid` output from the production image (wgrib2
|
||||||
|
/// 3.8.0) for an MSC HRDPS file. This is the shape the pipeline
|
||||||
|
/// actually sees; parsing only the `grid_template=10` form made every
|
||||||
|
/// HRDPS chain step fail with "could not parse nx×ny", which emptied
|
||||||
|
/// the Canadian half of the weather overlay.
|
||||||
|
#[test]
|
||||||
|
fn parse_wgrib2_grid_extracts_all_parameters_from_template_1() {
|
||||||
|
let output = "\
|
||||||
|
1:0:grid_template=1:winds(grid):
|
||||||
|
\trotated lat-lon grid:(2540 x 1290) units 1e-06 input WE:SN output WE:SN res 56
|
||||||
|
\tlat -12.302501 to 16.700001 by 0.022500
|
||||||
|
\tlon 345.178780 to 42.306283 by 0.022500 #points=3276600
|
||||||
|
\tsouth pole lat=-36.088520 lon=245.305142 angle of rot=0.000000
|
||||||
|
2:3593025:grid_template=1:winds(grid):
|
||||||
|
\trotated lat-lon grid:(2540 x 1290) units 1e-06 input WE:SN output WE:SN res 56
|
||||||
|
\tlat -12.302501 to 16.700001 by 0.022500
|
||||||
|
\tlon 345.178780 to 42.306283 by 0.022500 #points=3276600
|
||||||
|
\tsouth pole lat=-36.088520 lon=245.305142 angle of rot=0.000000
|
||||||
|
";
|
||||||
|
let params = parse_wgrib2_grid(output).unwrap();
|
||||||
|
assert_eq!(params.nx, 2540);
|
||||||
|
assert_eq!(params.ny, 1290);
|
||||||
|
assert!((params.sp_lat - -36.088520).abs() < 1e-6);
|
||||||
|
assert!((params.sp_lon - 245.305142).abs() < 1e-6);
|
||||||
|
assert!((params.lat0 - -12.302501).abs() < 1e-6);
|
||||||
|
assert!((params.lon0 - 345.178780).abs() < 1e-6);
|
||||||
|
assert!((params.d - 0.0225).abs() < 1e-6);
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The two templates describe the same physical grid, so they must
|
||||||
|
/// parse to (numerically) the same parameters.
|
||||||
|
#[test]
|
||||||
|
fn both_templates_agree_on_the_hrdps_grid() {
|
||||||
|
let template_1 = "\
|
||||||
|
1:0:grid_template=1:winds(grid):
|
||||||
|
\trotated lat-lon grid:(2540 x 1290) units 1e-06 input WE:SN output WE:SN res 56
|
||||||
|
\tlat -12.302501 to 16.700001 by 0.022500
|
||||||
|
\tlon 345.178780 to 42.306283 by 0.022500 #points=3276600
|
||||||
|
\tsouth pole lat=-36.088520 lon=245.305142 angle of rot=0.000000
|
||||||
|
";
|
||||||
|
let parsed = parse_wgrib2_grid(template_1).unwrap();
|
||||||
|
let expected = known_params();
|
||||||
|
|
||||||
|
assert_eq!(parsed.nx, expected.nx);
|
||||||
|
assert_eq!(parsed.ny, expected.ny);
|
||||||
|
assert!((parsed.sp_lat - expected.sp_lat).abs() < 1e-5);
|
||||||
|
assert!((parsed.sp_lon - expected.sp_lon).abs() < 1e-5);
|
||||||
|
assert!((parsed.lat0 - expected.lat0).abs() < 1e-5);
|
||||||
|
assert!((parsed.lon0 - expected.lon0).abs() < 1e-5);
|
||||||
|
assert!((parsed.d - expected.d).abs() < 1e-9);
|
||||||
|
}
|
||||||
|
|
||||||
|
#[test]
|
||||||
|
fn parse_wgrib2_grid_rejects_an_unrecognised_format() {
|
||||||
|
let output = "1:0:grid_template=0:winds(N/S):\n\tlat-lon grid:(1799 x 1059)\n";
|
||||||
|
assert!(parse_wgrib2_grid(output).is_err());
|
||||||
|
}
|
||||||
|
|
||||||
|
/// `RotatedPoleParams` has a single isotropic step; an anisotropic
|
||||||
|
/// grid must fail loudly rather than index rows with the column step.
|
||||||
|
#[test]
|
||||||
|
fn parse_wgrib2_grid_rejects_anisotropic_template_1_grid() {
|
||||||
|
let output = "\
|
||||||
|
1:0:grid_template=1:winds(grid):
|
||||||
|
\trotated lat-lon grid:(2540 x 1290) units 1e-06 input WE:SN output WE:SN res 56
|
||||||
|
\tlat -12.302501 to 16.700001 by 0.045000
|
||||||
|
\tlon 345.178780 to 42.306283 by 0.022500 #points=3276600
|
||||||
|
\tsouth pole lat=-36.088520 lon=245.305142 angle of rot=0.000000
|
||||||
|
";
|
||||||
|
assert!(parse_wgrib2_grid(output).is_err());
|
||||||
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
fn geo_to_rotated_returns_indices_for_known_point() {
|
fn geo_to_rotated_returns_indices_for_known_point() {
|
||||||
let params = known_params();
|
let params = known_params();
|
||||||
|
|
|
||||||
Loading…
Add table
Reference in a new issue