diff --git a/rust/prop_grid_rs/src/rotated_pole.rs b/rust/prop_grid_rs/src/rotated_pole.rs index 55853a9d..721da680 100644 --- a/rust/prop_grid_rs/src/rotated_pole.rs +++ b/rust/prop_grid_rs/src/rotated_pole.rs @@ -1,8 +1,9 @@ //! Rotated lat/lon → geographic coordinate transform for HRDPS. //! -//! HRDPS uses a rotated-pole grid (GRIB2 template 10) where the -//! computational grid is regular lat/lon but the geographic grid is -//! rotated. `wgrib2 -lola` handles this by brute-forcing per-output-point, +//! HRDPS uses a rotated-pole grid (GRIB2 grid template 1; wgrib2 also +//! renders the equivalent template 10) where the computational grid is +//! 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 //! bbox at 0.125°. //! @@ -175,7 +176,19 @@ fn build_lookup(params: &RotatedPoleParams, target_spec: &GridSpec) -> Vec Vec Result { - // Extract nx, ny from "Rotated Lat/lon Grid: (NNNN x NNNN)" - let (nx, ny) = output - .lines() - .find_map(|line| { - 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 - } + parse_template_10(output) + .or_else(|| parse_template_1(output)) + .ok_or_else(|| { + format!("could not parse rotated-pole grid from wgrib2 -grid output:\n{output}") }) - .ok_or_else(|| format!("could not parse nx×ny from wgrib2 -grid output:\n{output}"))?; +} - // Extract lat0 (LatFirst), lon0 (LonFirst), d (Di) - let lat0 = parse_float_after(output, "LatFirst:").ok_or("LatFirst not found")?; - let lon0 = parse_float_after(output, "LonFirst:").ok_or("LonFirst not found")?; - let d = parse_float_after(output, "Di:").ok_or("Di not found")?; +/// Verbose `grid_template=10` form: one labelled `key: value` per line. +fn parse_template_10(output: &str) -> Option { + let (nx, ny) = output.lines().find_map(|line| { + let rest = line.trim().strip_prefix("Rotated Lat/lon Grid:")?; + parse_dims(rest.trim()) + })?; - // Extract south pole location - let sp_lat = output + let lat0 = parse_float_after(output, "LatFirst:")?; + let lon0 = parse_float_after(output, "LonFirst:")?; + let d = parse_float_after(output, "Di:")?; + + let pole_line = output .lines() .skip_while(|l| !l.trim().starts_with("South Pole Location:")) - .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")?; + .nth(1)?; - Ok(RotatedPoleParams { - sp_lat, - sp_lon, + Some(RotatedPoleParams { + sp_lat: parse_float_after(pole_line, "lat:")?, + sp_lon: parse_float_after(pole_line, "lon:")?, lat0, lon0, d, @@ -239,6 +238,64 @@ fn parse_wgrib2_grid(output: &str) -> Result { }) } +/// Compact `grid_template=1` form. Each axis line carries the origin and +/// the step together (`lat to by `), and the pole +/// sits on a single `key=value` line. +fn parse_template_1(output: &str) -> Option { + 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 ` to by ` 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 { for line in haystack.lines() { if let Some(pos) = line.find(prefix) { @@ -290,6 +347,78 @@ mod tests { 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] fn geo_to_rotated_returns_indices_for_known_point() { let params = known_params();