From e3be41e6c2bf8cd047d124e89547064ac7c9afe8 Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 04:51:56 -0500 Subject: [PATCH 1/7] Decode rotated and regular Gaussian grids and reject quasi-regular geometry --- CHANGELOG.md | 2 + README.md | 5 +- grib-core/src/error.rs | 12 + grib-core/src/grid.rs | 704 +++++++++++++++++++++++++-- grib-core/src/lib.rs | 2 +- grib-reader/src/grid.rs | 2 +- grib-reader/src/lib.rs | 129 ++++- grib-reader/tests/common/fixtures.rs | 62 +++ grib-reader/tests/common/mod.rs | 3 +- grib-reader/tests/parity_eccodes.rs | 11 + 10 files changed, 890 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 670cedf..87f3eb6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,8 @@ expand bitmap gaps a byte at a time - locate GRIB messages with substring search instead of testing every byte offset +- decode rotated latitude/longitude and regular Gaussian grids while rejecting + quasi-regular geometry with a typed error - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index 75eb882..0e12192 100644 --- a/README.md +++ b/README.md @@ -163,6 +163,8 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Logical field indexing for multi-field GRIB2 messages - GRIB2 Section 2 access and same-message bitmap reuse through indicator 254 - Regular latitude/longitude grids for GRIB1 and GRIB2 +- Reader GRIB2 rotated latitude/longitude template 3.1 and regular Gaussian + template 3.40, with explicit typed rejection of quasi-regular grids - Reader GRIB2 Mercator grid template 3.10, polar stereographic grid template 3.20, Lambert conformal grid template 3.30, and Albers equal-area grid template 3.31 metadata, projected coordinate offsets, and flat data decode @@ -194,7 +196,8 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; ## Not Yet Supported -- Remaining GRIB2 grid templates beyond 3.0, 3.10, 3.20, 3.30, and 3.31 +- Remaining GRIB2 grid templates beyond 3.0, 3.1, 3.10, 3.20, 3.30, 3.31, + and 3.40 - Remaining GRIB2 product definition templates beyond 4.0, 4.1, 4.8, and 4.11 - Writer GRIB2 row-by-row complex packing diff --git a/grib-core/src/error.rs b/grib-core/src/error.rs index 8827de5..65df93d 100644 --- a/grib-core/src/error.rs +++ b/grib-core/src/error.rs @@ -29,6 +29,18 @@ pub enum Error { #[error("unsupported grid definition template: {0}")] UnsupportedGridTemplate(u16), + #[error("unsupported quasi-regular grid definition template: {template}")] + UnsupportedQuasiRegularGrid { template: u16 }, + + #[error( + "grid template {template} uses unsupported angular unit {basic_angle}/{subdivisions} degrees" + )] + UnsupportedGridAngularUnit { + template: u16, + basic_angle: u32, + subdivisions: u32, + }, + #[error("unsupported data representation template: {0}")] UnsupportedDataTemplate(u16), diff --git a/grib-core/src/grid.rs b/grib-core/src/grid.rs index dc4e03c..eecbf3c 100644 --- a/grib-core/src/grid.rs +++ b/grib-core/src/grid.rs @@ -1,8 +1,8 @@ //! Grid Definition Section (Section 3) parsing. use crate::binary::decode_wmo_i32; -use crate::ensure_limit; use crate::error::{Error, Result}; +use crate::{ensure_limit, filled_vec}; /// Grid definition extracted from Section 3. /// @@ -15,6 +15,8 @@ use crate::error::{Error, Result}; pub enum GridDefinition { /// Template 3.0: Regular latitude/longitude (equidistant cylindrical). LatLon(LatLonGrid), + /// Template 3.1: Rotated latitude/longitude. + RotatedLatLon(RotatedLatLonGrid), /// Template 3.10: Mercator. Mercator(MercatorGrid), /// Template 3.31: Albers equal-area. @@ -23,6 +25,8 @@ pub enum GridDefinition { LambertConformal(LambertConformalGrid), /// Template 3.20: Polar stereographic projection. PolarStereographic(PolarStereographicGrid), + /// Template 3.40: Regular Gaussian latitude/longitude. + RegularGaussian(RegularGaussianGrid), /// Unsupported template (stored for diagnostics). Unsupported(u16), } @@ -41,6 +45,31 @@ pub struct LatLonGrid { pub scanning_mode: u8, } +/// Template 3.1: rotated latitude/longitude grid. +#[derive(Debug, Clone, PartialEq)] +pub struct RotatedLatLonGrid { + /// Grid coordinates expressed in the rotated coordinate system. + pub grid: LatLonGrid, + pub lat_southern_pole: i32, + pub lon_southern_pole: u32, + pub angle_of_rotation: f32, +} + +/// Template 3.40: regular Gaussian latitude/longitude grid. +#[derive(Debug, Clone, PartialEq)] +pub struct RegularGaussianGrid { + pub ni: u32, + pub nj: u32, + pub lat_first: i32, + pub lon_first: i32, + pub lat_last: i32, + pub lon_last: i32, + pub di: u32, + /// Number of Gaussian parallels between a pole and the equator. + pub number_of_parallels: u32, + pub scanning_mode: u8, +} + /// Template 3.10: Mercator grid. #[derive(Debug, Clone, PartialEq)] pub struct MercatorGrid { @@ -119,10 +148,12 @@ impl GridDefinition { pub fn template_number(&self) -> u16 { match self { Self::LatLon(_) => 0, + Self::RotatedLatLon(_) => 1, Self::Mercator(_) => 10, Self::PolarStereographic(_) => 20, Self::LambertConformal(_) => 30, Self::AlbersEqualArea(_) => 31, + Self::RegularGaussian(_) => 40, Self::Unsupported(template) => *template, } } @@ -134,6 +165,20 @@ impl GridDefinition { } } + pub fn as_rotated_lat_lon(&self) -> Option<&RotatedLatLonGrid> { + match self { + Self::RotatedLatLon(grid) => Some(grid), + _ => None, + } + } + + pub fn as_regular_gaussian(&self) -> Option<&RegularGaussianGrid> { + match self { + Self::RegularGaussian(grid) => Some(grid), + _ => None, + } + } + pub fn as_polar_stereographic(&self) -> Option<&PolarStereographicGrid> { match self { Self::PolarStereographic(grid) => Some(grid), @@ -168,7 +213,10 @@ impl GridDefinition { Self::PolarStereographic(grid) => Some(&grid.core), Self::LambertConformal(grid) => Some(&grid.core), Self::AlbersEqualArea(grid) => Some(&grid.core), - Self::LatLon(_) | Self::Unsupported(_) => None, + Self::LatLon(_) + | Self::RotatedLatLon(_) + | Self::RegularGaussian(_) + | Self::Unsupported(_) => None, } } @@ -182,10 +230,12 @@ impl GridDefinition { pub fn shape(&self) -> (usize, usize) { match self { Self::LatLon(g) => (g.ni as usize, g.nj as usize), + Self::RotatedLatLon(g) => (g.grid.ni as usize, g.grid.nj as usize), Self::Mercator(g) => (g.core.nx as usize, g.core.ny as usize), Self::PolarStereographic(g) => (g.core.nx as usize, g.core.ny as usize), Self::LambertConformal(g) => (g.core.nx as usize, g.core.ny as usize), Self::AlbersEqualArea(g) => (g.core.nx as usize, g.core.ny as usize), + Self::RegularGaussian(g) => (g.ni as usize, g.nj as usize), Self::Unsupported(_) => (0, 0), } } @@ -193,10 +243,12 @@ impl GridDefinition { pub fn shape_num_points(&self) -> Result { match self { Self::LatLon(g) => checked_grid_point_count(g.ni, g.nj), + Self::RotatedLatLon(g) => checked_grid_point_count(g.grid.ni, g.grid.nj), Self::Mercator(g) => checked_grid_point_count(g.core.nx, g.core.ny), Self::PolarStereographic(g) => checked_grid_point_count(g.core.nx, g.core.ny), Self::LambertConformal(g) => checked_grid_point_count(g.core.nx, g.core.ny), Self::AlbersEqualArea(g) => checked_grid_point_count(g.core.nx, g.core.ny), + Self::RegularGaussian(g) => checked_grid_point_count(g.ni, g.nj), Self::Unsupported(_) => Ok(0), } } @@ -205,6 +257,8 @@ impl GridDefinition { let (ni, nj) = self.shape(); match self { Self::LatLon(_) + | Self::RotatedLatLon(_) + | Self::RegularGaussian(_) | Self::Mercator(_) | Self::PolarStereographic(_) | Self::LambertConformal(_) @@ -224,10 +278,12 @@ impl GridDefinition { pub fn checked_num_points(&self) -> Result { match self { Self::LatLon(_) => self.shape_num_points(), + Self::RotatedLatLon(_) => self.shape_num_points(), Self::Mercator(g) => Ok(g.core.number_of_points as usize), Self::PolarStereographic(g) => Ok(g.core.number_of_points as usize), Self::LambertConformal(g) => Ok(g.core.number_of_points as usize), Self::AlbersEqualArea(g) => Ok(g.core.number_of_points as usize), + Self::RegularGaussian(_) => self.shape_num_points(), Self::Unsupported(_) => Ok(0), } } @@ -238,7 +294,39 @@ impl GridDefinition { Self::PolarStereographic(g) => Some(g.core.number_of_points as usize), Self::LambertConformal(g) => Some(g.core.number_of_points as usize), Self::AlbersEqualArea(g) => Some(g.core.number_of_points as usize), - Self::LatLon(_) | Self::Unsupported(_) => None, + Self::LatLon(_) + | Self::RotatedLatLon(_) + | Self::RegularGaussian(_) + | Self::Unsupported(_) => None, + } + } + + /// One-dimensional latitude axis for unrotated geographic grids. + pub fn latitudes(&self) -> Result>> { + self.latitudes_with_limit(None) + } + + pub fn latitudes_with_limit(&self, max_axis_points: Option) -> Result>> { + match self { + Self::LatLon(grid) => Ok(Some(grid.latitudes_with_limit(max_axis_points)?)), + Self::RegularGaussian(grid) => Ok(Some(grid.latitudes_with_limit(max_axis_points)?)), + _ => Ok(None), + } + } + + /// One-dimensional longitude axis for unrotated geographic grids. + pub fn longitudes(&self) -> Result>> { + self.longitudes_with_limit(None) + } + + pub fn longitudes_with_limit( + &self, + max_axis_points: Option, + ) -> Result>> { + match self { + Self::LatLon(grid) => Ok(Some(grid.longitudes_with_limit(max_axis_points)?)), + Self::RegularGaussian(grid) => Ok(Some(grid.longitudes_with_limit(max_axis_points)?)), + _ => Ok(None), } } @@ -299,10 +387,12 @@ impl GridDefinition { pub fn validate_supported_scan_order(&self) -> Result<()> { match self { Self::LatLon(grid) => grid.validate_supported_scan_order(), + Self::RotatedLatLon(grid) => grid.grid.validate_supported_scan_order(), Self::Mercator(grid) => grid.core.validate_supported_scan_order(), Self::PolarStereographic(grid) => grid.core.validate_supported_scan_order(), Self::LambertConformal(grid) => grid.core.validate_supported_scan_order(), Self::AlbersEqualArea(grid) => grid.core.validate_supported_scan_order(), + Self::RegularGaussian(grid) => grid.validate_supported_scan_order(), Self::Unsupported(template) => Err(Error::UnsupportedGridTemplate(*template)), } } @@ -312,10 +402,17 @@ impl GridDefinition { pub fn reorder_for_ndarray_in_place(&self, values: &mut [T]) -> Result<()> { match self { Self::LatLon(grid) => grid.reorder_for_ndarray_in_place(values), + Self::RotatedLatLon(grid) => grid.grid.reorder_for_ndarray_in_place(values), Self::Mercator(grid) => grid.core.reorder_for_ndarray_in_place(values), Self::PolarStereographic(grid) => grid.core.reorder_for_ndarray_in_place(values), Self::LambertConformal(grid) => grid.core.reorder_for_ndarray_in_place(values), Self::AlbersEqualArea(grid) => grid.core.reorder_for_ndarray_in_place(values), + Self::RegularGaussian(grid) => transform_supported_scan_order_in_place( + values, + grid.ni as usize, + grid.nj as usize, + grid.scanning_mode, + ), Self::Unsupported(template) => Err(Error::UnsupportedGridTemplate(*template)), } } @@ -331,6 +428,12 @@ impl GridDefinition { grid.nj as usize, grid.scanning_mode, ), + Self::RotatedLatLon(grid) => normalize_north_up_in_place( + values, + grid.grid.ni as usize, + grid.grid.nj as usize, + grid.grid.scanning_mode, + ), Self::Mercator(grid) => normalize_north_up_in_place( values, grid.core.nx as usize, @@ -355,6 +458,12 @@ impl GridDefinition { grid.core.ny as usize, grid.core.scanning_mode, ), + Self::RegularGaussian(grid) => normalize_north_up_in_place( + values, + grid.ni as usize, + grid.nj as usize, + grid.scanning_mode, + ), Self::Unsupported(template) => Err(Error::UnsupportedGridTemplate(*template)), } } @@ -376,10 +485,12 @@ impl GridDefinition { let template = u16::from_be_bytes(section_bytes[12..14].try_into().unwrap()); match template { 0 => parse_latlon(section_bytes), + 1 => parse_rotated_latlon(section_bytes), 10 => parse_mercator(section_bytes), 20 => parse_polar_stereographic(section_bytes), 30 => parse_lambert_conformal(section_bytes), 31 => parse_albers_equal_area(section_bytes), + 40 => parse_regular_gaussian(section_bytes), _ => Ok(Self::Unsupported(template)), } } @@ -470,6 +581,236 @@ impl LatLonGrid { } } +impl RotatedLatLonGrid { + pub fn rotated_longitudes(&self) -> Result> { + self.grid.longitudes() + } + + pub fn rotated_longitudes_with_limit( + &self, + max_axis_points: Option, + ) -> Result> { + self.grid.longitudes_with_limit(max_axis_points) + } + + pub fn rotated_latitudes(&self) -> Result> { + self.grid.latitudes() + } + + pub fn rotated_latitudes_with_limit(&self, max_axis_points: Option) -> Result> { + self.grid.latitudes_with_limit(max_axis_points) + } +} + +impl RegularGaussianGrid { + pub fn longitudes(&self) -> Result> { + self.longitudes_with_limit(None) + } + + pub fn longitudes_with_limit(&self, max_axis_points: Option) -> Result> { + let step = self.longitude_increment_degrees()?; + let signed_step = if i_scans_positive(self.scanning_mode) { + step + } else { + -step + }; + linear_axis( + self.ni, + f64::from(self.lon_first) / 1_000_000.0, + signed_step, + max_axis_points, + "Gaussian longitude axis", + ) + } + + pub fn latitudes(&self) -> Result> { + self.latitudes_with_limit(None) + } + + pub fn latitudes_with_limit(&self, max_axis_points: Option) -> Result> { + let requested = self.nj as usize; + ensure_limit("Gaussian latitude axis", requested, max_axis_points)?; + let all_latitudes = gaussian_latitudes(self.number_of_parallels, max_axis_points)?; + if requested > all_latitudes.len() { + return Err(Error::InvalidSection { + section: 3, + reason: format!( + "Gaussian grid requests {} rows but N={} defines only {} parallels", + self.nj, + self.number_of_parallels, + all_latitudes.len() + ), + }); + } + if requested == 0 { + return Ok(Vec::new()); + } + + let first = f64::from(self.lat_first) / 1_000_000.0; + let first_index = all_latitudes + .iter() + .enumerate() + .min_by(|(_, left), (_, right)| { + (*left - first).abs().total_cmp(&(*right - first).abs()) + }) + .map(|(index, _)| index) + .ok_or_else(|| Error::InvalidSection { + section: 3, + reason: "Gaussian grid has no latitude parallels".into(), + })?; + if (all_latitudes[first_index] - first).abs() > 0.001 { + return Err(Error::InvalidSection { + section: 3, + reason: format!( + "first latitude {first} is not a Gaussian parallel for N={}", + self.number_of_parallels + ), + }); + } + + let ascending = j_scans_positive(self.scanning_mode); + let last_index = if ascending { + first_index.checked_sub(requested - 1) + } else { + first_index.checked_add(requested - 1) + } + .filter(|index| *index < all_latitudes.len()) + .ok_or_else(|| Error::InvalidSection { + section: 3, + reason: "Gaussian latitude scan extends beyond the defined parallels".into(), + })?; + + let mut result = Vec::new(); + result.try_reserve_exact(requested).map_err(|error| { + Error::allocation("Gaussian latitude coordinates", requested, error) + })?; + if ascending { + result.extend( + (last_index..=first_index) + .rev() + .map(|index| all_latitudes[index]), + ); + } else { + result.extend(all_latitudes[first_index..=last_index].iter().copied()); + } + + let declared_last = f64::from(self.lat_last) / 1_000_000.0; + if (result[requested - 1] - declared_last).abs() > 0.001 { + return Err(Error::InvalidSection { + section: 3, + reason: format!( + "last latitude {declared_last} does not match the Gaussian scan endpoint {}", + result[requested - 1] + ), + }); + } + Ok(result) + } + + pub fn validate_supported_scan_order(&self) -> Result<()> { + validate_supported_scan_order(self.scanning_mode) + } + + fn longitude_increment_degrees(&self) -> Result { + if self.di != u32::MAX { + let increment = f64::from(self.di) / 1_000_000.0; + self.validate_longitude_endpoint(increment)?; + return Ok(increment); + } + if self.ni == 0 { + return Err(Error::InvalidSection { + section: 3, + reason: "Gaussian grid has missing Di and zero Ni".into(), + }); + } + let inferred = 360.0 / f64::from(self.ni); + self.validate_longitude_endpoint(inferred)?; + Ok(inferred) + } + + fn validate_longitude_endpoint(&self, increment: f64) -> Result<()> { + if self.ni == 0 { + return Err(Error::InvalidSection { + section: 3, + reason: "Gaussian grid has zero Ni".into(), + }); + } + let expected_last = f64::from(self.lon_first) / 1_000_000.0 + + if i_scans_positive(self.scanning_mode) { + increment * f64::from(self.ni - 1) + } else { + -increment * f64::from(self.ni - 1) + }; + let declared_last = f64::from(self.lon_last) / 1_000_000.0; + let wrapped_difference = (expected_last - declared_last + 180.0).rem_euclid(360.0) - 180.0; + if wrapped_difference.abs() > 0.001 { + return Err(Error::InvalidSection { + section: 3, + reason: "Gaussian longitude endpoint does not match Ni, Di, and scanning mode" + .into(), + }); + } + Ok(()) + } +} + +fn gaussian_latitudes( + number_of_parallels: u32, + max_axis_points: Option, +) -> Result> { + if number_of_parallels == 0 { + return Err(Error::InvalidSection { + section: 3, + reason: "Gaussian grid N must be nonzero".into(), + }); + } + let order = usize::try_from(u64::from(number_of_parallels) * 2) + .map_err(|_| Error::ValueOutOfRange("Gaussian parallel count does not fit usize".into()))?; + ensure_limit("Gaussian parallel calculation", order, max_axis_points)?; + let mut latitudes = filled_vec(order, 0.0, "Gaussian latitude coordinates")?; + let order_f64 = order as f64; + + for index in 0..order / 2 { + let mut root = (std::f64::consts::PI * (index as f64 + 0.75) / (order_f64 + 0.5)).cos(); + let mut converged = false; + for _ in 0..32 { + let (polynomial, previous) = legendre_pair(order, root); + let derivative = order_f64 * (root * polynomial - previous) / (root * root - 1.0); + let next = root - polynomial / derivative; + if (next - root).abs() <= 4.0 * f64::EPSILON * next.abs().max(1.0) { + root = next; + converged = true; + break; + } + root = next; + } + if !converged || !root.is_finite() { + return Err(Error::InvalidSection { + section: 3, + reason: format!( + "failed to compute Gaussian latitude root {index} for N={number_of_parallels}" + ), + }); + } + let latitude = root.asin().to_degrees(); + latitudes[index] = latitude; + latitudes[order - index - 1] = -latitude; + } + Ok(latitudes) +} + +fn legendre_pair(order: usize, value: f64) -> (f64, f64) { + let mut previous = 0.0; + let mut current = 1.0; + for degree in 1..=order { + let older = previous; + previous = current; + current = ((2 * degree - 1) as f64 * value * previous - (degree - 1) as f64 * older) + / degree as f64; + } + (current, previous) +} + impl ProjectedGridCore { pub fn x_coordinates(&self) -> Result> { self.x_coordinates_with_limit(None) @@ -508,9 +849,9 @@ fn transform_supported_scan_order_in_place( scanning_mode: u8, ) -> Result<()> { validate_supported_scan_order(scanning_mode)?; - let expected = ni - .checked_mul(nj) - .ok_or_else(|| Error::Other("grid point count overflow".into()))?; + let expected = ni.checked_mul(nj).ok_or(Error::ArithmeticOverflow { + operation: "computing grid point count", + })?; if values.len() != expected { return Err(Error::DataLengthMismatch { expected, @@ -637,8 +978,9 @@ fn linear_axis( fn checked_grid_point_count(nx: u32, ny: u32) -> Result { let count = u64::from(nx) * u64::from(ny); - usize::try_from(count) - .map_err(|_| Error::Other(format!("grid point count {count} does not fit in usize"))) + usize::try_from(count).map_err(|_| { + Error::ValueOutOfRange(format!("grid point count {count} does not fit in usize")) + }) } fn reverse_alternating_rows(values: &mut [T], ni: usize, nj: usize, i_scans_positive: bool) { @@ -662,27 +1004,171 @@ fn parse_latlon(data: &[u8]) -> Result { }); } + reject_quasi_regular_grid(data, 0)?; + require_microdegree_angular_unit(data, 0)?; + let grid = parse_latlon_fields(data); + validate_geographic_coordinates( + grid.lat_first, + grid.lon_first, + grid.lat_last, + grid.lon_last, + 0, + )?; + Ok(GridDefinition::LatLon(grid)) +} + +fn parse_rotated_latlon(data: &[u8]) -> Result { + if data.len() < 84 { + return Err(Error::InvalidSection { + section: 3, + reason: format!("template 3.1 requires 84 bytes, got {}", data.len()), + }); + } + + reject_quasi_regular_grid(data, 1)?; + require_microdegree_angular_unit(data, 1)?; + let grid = parse_latlon_fields(data); + validate_geographic_coordinates( + grid.lat_first, + grid.lon_first, + grid.lat_last, + grid.lon_last, + 1, + )?; + + let lat_southern_pole = decode_wmo_i32(&data[72..76]).unwrap(); + let lon_southern_pole = u32::from_be_bytes(data[76..80].try_into().unwrap()); + if lat_southern_pole.unsigned_abs() > 90_000_000 || lon_southern_pole > 360_000_000 { + return Err(Error::InvalidSection { + section: 3, + reason: "rotated-grid southern pole is outside geographic bounds".into(), + }); + } + let angle_of_rotation = f32::from_be_bytes(data[80..84].try_into().unwrap()); + if !angle_of_rotation.is_finite() { + return Err(Error::InvalidSection { + section: 3, + reason: "rotated-grid angle of rotation must be finite".into(), + }); + } + + Ok(GridDefinition::RotatedLatLon(RotatedLatLonGrid { + grid, + lat_southern_pole, + lon_southern_pole, + angle_of_rotation, + })) +} + +fn parse_regular_gaussian(data: &[u8]) -> Result { + if data.len() < 72 { + return Err(Error::InvalidSection { + section: 3, + reason: format!("template 3.40 requires 72 bytes, got {}", data.len()), + }); + } + + reject_quasi_regular_grid(data, 40)?; + require_microdegree_angular_unit(data, 40)?; + let grid = RegularGaussianGrid { + ni: u32::from_be_bytes(data[30..34].try_into().unwrap()), + nj: u32::from_be_bytes(data[34..38].try_into().unwrap()), + lat_first: decode_wmo_i32(&data[46..50]).unwrap(), + lon_first: decode_wmo_i32(&data[50..54]).unwrap(), + lat_last: decode_wmo_i32(&data[55..59]).unwrap(), + lon_last: decode_wmo_i32(&data[59..63]).unwrap(), + di: u32::from_be_bytes(data[63..67].try_into().unwrap()), + number_of_parallels: u32::from_be_bytes(data[67..71].try_into().unwrap()), + scanning_mode: data[71], + }; + validate_geographic_coordinates( + grid.lat_first, + grid.lon_first, + grid.lat_last, + grid.lon_last, + 40, + )?; + if grid.number_of_parallels == 0 { + return Err(Error::InvalidSection { + section: 3, + reason: "template 3.40 number of parallels N must be nonzero".into(), + }); + } + let available_rows = u64::from(grid.number_of_parallels) * 2; + if u64::from(grid.nj) > available_rows { + return Err(Error::InvalidSection { + section: 3, + reason: format!( + "template 3.40 Nj={} exceeds the {available_rows} parallels defined by N={}", + grid.nj, grid.number_of_parallels + ), + }); + } + Ok(GridDefinition::RegularGaussian(grid)) +} + +fn parse_latlon_fields(data: &[u8]) -> LatLonGrid { + LatLonGrid { + ni: u32::from_be_bytes(data[30..34].try_into().unwrap()), + nj: u32::from_be_bytes(data[34..38].try_into().unwrap()), + lat_first: decode_wmo_i32(&data[46..50]).unwrap(), + lon_first: decode_wmo_i32(&data[50..54]).unwrap(), + lat_last: decode_wmo_i32(&data[55..59]).unwrap(), + lon_last: decode_wmo_i32(&data[59..63]).unwrap(), + di: u32::from_be_bytes(data[63..67].try_into().unwrap()), + dj: u32::from_be_bytes(data[67..71].try_into().unwrap()), + scanning_mode: data[71], + } +} + +fn reject_quasi_regular_grid(data: &[u8], template: u16) -> Result<()> { let ni = u32::from_be_bytes(data[30..34].try_into().unwrap()); let nj = u32::from_be_bytes(data[34..38].try_into().unwrap()); - let lat_first = decode_wmo_i32(&data[46..50]).unwrap(); - let lon_first = decode_wmo_i32(&data[50..54]).unwrap(); - let lat_last = decode_wmo_i32(&data[55..59]).unwrap(); - let lon_last = decode_wmo_i32(&data[59..63]).unwrap(); - let di = u32::from_be_bytes(data[63..67].try_into().unwrap()); - let dj = u32::from_be_bytes(data[67..71].try_into().unwrap()); - let scanning_mode = data[71]; - - Ok(GridDefinition::LatLon(LatLonGrid { - ni, - nj, - lat_first, - lon_first, - lat_last, - lon_last, - di, - dj, - scanning_mode, - })) + if ni == u32::MAX || nj == u32::MAX || data[10] != 0 { + return Err(Error::UnsupportedQuasiRegularGrid { template }); + } + Ok(()) +} + +fn require_microdegree_angular_unit(data: &[u8], template: u16) -> Result<()> { + let encoded_basic_angle = u32::from_be_bytes(data[38..42].try_into().unwrap()); + let encoded_subdivisions = u32::from_be_bytes(data[42..46].try_into().unwrap()); + let basic_angle = match encoded_basic_angle { + 0 | u32::MAX => 1, + value => value, + }; + let subdivisions = match encoded_subdivisions { + 0 | u32::MAX => 1_000_000, + value => value, + }; + if u64::from(basic_angle) * 1_000_000 != u64::from(subdivisions) { + return Err(Error::UnsupportedGridAngularUnit { + template, + basic_angle, + subdivisions, + }); + } + Ok(()) +} + +fn validate_geographic_coordinates( + lat_first: i32, + lon_first: i32, + lat_last: i32, + lon_last: i32, + template: u16, +) -> Result<()> { + if lat_first.unsigned_abs() > 90_000_000 + || lat_last.unsigned_abs() > 90_000_000 + || lon_first.unsigned_abs() > 360_000_000 + || lon_last.unsigned_abs() > 360_000_000 + { + return Err(Error::InvalidSection { + section: 3, + reason: format!("template 3.{template} coordinates are outside geographic bounds"), + }); + } + Ok(()) } fn parse_mercator(data: &[u8]) -> Result { @@ -791,7 +1277,7 @@ fn parse_projected_core( mod tests { use super::{ AlbersEqualAreaGrid, GridDefinition, LambertConformalGrid, LatLonGrid, MercatorGrid, - PolarStereographicGrid, ProjectedGridCore, + PolarStereographicGrid, ProjectedGridCore, RegularGaussianGrid, }; use crate::binary::encode_wmo_i32; @@ -882,6 +1368,130 @@ mod tests { } } + #[test] + fn parses_rotated_latlon_template() { + let section = build_rotated_latlon_section(); + let grid = GridDefinition::parse(§ion).unwrap(); + + assert_eq!(grid.shape(), (3, 2)); + assert_eq!(grid.ndarray_shape(), vec![2, 3]); + assert_eq!(grid.template_number(), 1); + assert!(grid.as_lat_lon().is_none()); + let rotated = grid.as_rotated_lat_lon().unwrap(); + assert_eq!(rotated.lat_southern_pole, -30_000_000); + assert_eq!(rotated.lon_southern_pole, 10_000_000); + assert_eq!(rotated.angle_of_rotation, 15.5); + assert_eq!( + rotated.rotated_longitudes().unwrap(), + vec![-20.0, 0.0, 20.0] + ); + assert_eq!(rotated.rotated_latitudes().unwrap(), vec![-10.0, 0.0]); + assert_eq!(grid.latitudes().unwrap(), None); + } + + #[test] + fn rejects_non_finite_rotated_grid_angle() { + let mut section = build_rotated_latlon_section(); + section[80..84].copy_from_slice(&f32::NAN.to_be_bytes()); + assert!(matches!( + GridDefinition::parse(§ion), + Err(crate::Error::InvalidSection { section: 3, .. }) + )); + } + + #[test] + fn parses_regular_gaussian_and_computes_latitudes() { + let section = build_regular_gaussian_section(90_000_000); + let grid = GridDefinition::parse(§ion).unwrap(); + + assert_eq!(grid.shape(), (4, 4)); + assert_eq!(grid.ndarray_shape(), vec![4, 4]); + assert_eq!(grid.template_number(), 40); + let gaussian = grid.as_regular_gaussian().unwrap(); + assert_eq!( + gaussian, + &RegularGaussianGrid { + ni: 4, + nj: 4, + lat_first: 59_444_408, + lon_first: 0, + lat_last: -59_444_408, + lon_last: 270_000_000, + di: 90_000_000, + number_of_parallels: 2, + scanning_mode: 0, + } + ); + let latitudes = gaussian.latitudes().unwrap(); + let expected = [ + 59.444_408_289_166_77, + 19.875_719_147_440_904, + -19.875_719_147_440_904, + -59.444_408_289_166_77, + ]; + for (actual, expected) in latitudes.iter().zip(expected) { + assert!((actual - expected).abs() < 1e-12); + } + assert_eq!( + gaussian.longitudes().unwrap(), + vec![0.0, 90.0, 180.0, 270.0] + ); + assert_eq!(grid.latitudes().unwrap().unwrap(), latitudes); + } + + #[test] + fn infers_missing_gaussian_increment_for_global_grid() { + let section = build_regular_gaussian_section(u32::MAX); + let GridDefinition::RegularGaussian(grid) = GridDefinition::parse(§ion).unwrap() else { + panic!("expected regular Gaussian grid"); + }; + assert_eq!(grid.longitudes().unwrap(), vec![0.0, 90.0, 180.0, 270.0]); + } + + #[test] + fn rejects_inconsistent_gaussian_longitude_endpoint() { + let mut section = build_regular_gaussian_section(90_000_000); + section[59..63].copy_from_slice(&encode_wmo_i32(180_000_000).unwrap()); + let GridDefinition::RegularGaussian(grid) = GridDefinition::parse(§ion).unwrap() else { + panic!("expected regular Gaussian grid"); + }; + assert!(matches!( + grid.longitudes(), + Err(crate::Error::InvalidSection { section: 3, .. }) + )); + } + + #[test] + fn rejects_quasi_regular_geographic_grids_with_typed_error() { + for (template, mut section) in [ + (1, build_rotated_latlon_section()), + (40, build_regular_gaussian_section(90_000_000)), + ] { + section[30..34].copy_from_slice(&u32::MAX.to_be_bytes()); + section[10] = 2; + assert!(matches!( + GridDefinition::parse(§ion), + Err(crate::Error::UnsupportedQuasiRegularGrid { template: actual }) + if actual == template + )); + } + } + + #[test] + fn rejects_non_microdegree_geographic_angle_units_explicitly() { + let mut section = build_regular_gaussian_section(90_000_000); + section[38..42].copy_from_slice(&1u32.to_be_bytes()); + section[42..46].copy_from_slice(&1_000u32.to_be_bytes()); + assert!(matches!( + GridDefinition::parse(§ion), + Err(crate::Error::UnsupportedGridAngularUnit { + template: 40, + basic_angle: 1, + subdivisions: 1_000, + }) + )); + } + #[test] fn parses_polar_stereographic_template() { let section = build_polar_stereographic_section(0); @@ -1279,6 +1889,44 @@ mod tests { section } + fn build_rotated_latlon_section() -> Vec { + let mut section = vec![0u8; 84]; + section[..4].copy_from_slice(&84u32.to_be_bytes()); + section[4] = 3; + section[6..10].copy_from_slice(&6u32.to_be_bytes()); + section[12..14].copy_from_slice(&1u16.to_be_bytes()); + section[30..34].copy_from_slice(&3u32.to_be_bytes()); + section[34..38].copy_from_slice(&2u32.to_be_bytes()); + section[46..50].copy_from_slice(&encode_wmo_i32(-10_000_000).unwrap()); + section[50..54].copy_from_slice(&encode_wmo_i32(-20_000_000).unwrap()); + section[55..59].copy_from_slice(&encode_wmo_i32(0).unwrap()); + section[59..63].copy_from_slice(&encode_wmo_i32(20_000_000).unwrap()); + section[63..67].copy_from_slice(&20_000_000u32.to_be_bytes()); + section[67..71].copy_from_slice(&10_000_000u32.to_be_bytes()); + section[71] = 0b0100_0000; + section[72..76].copy_from_slice(&encode_wmo_i32(-30_000_000).unwrap()); + section[76..80].copy_from_slice(&10_000_000u32.to_be_bytes()); + section[80..84].copy_from_slice(&15.5f32.to_be_bytes()); + section + } + + fn build_regular_gaussian_section(di: u32) -> Vec { + let mut section = vec![0u8; 72]; + section[..4].copy_from_slice(&72u32.to_be_bytes()); + section[4] = 3; + section[6..10].copy_from_slice(&16u32.to_be_bytes()); + section[12..14].copy_from_slice(&40u16.to_be_bytes()); + section[30..34].copy_from_slice(&4u32.to_be_bytes()); + section[34..38].copy_from_slice(&4u32.to_be_bytes()); + section[46..50].copy_from_slice(&encode_wmo_i32(59_444_408).unwrap()); + section[50..54].copy_from_slice(&encode_wmo_i32(0).unwrap()); + section[55..59].copy_from_slice(&encode_wmo_i32(-59_444_408).unwrap()); + section[59..63].copy_from_slice(&encode_wmo_i32(270_000_000).unwrap()); + section[63..67].copy_from_slice(&di.to_be_bytes()); + section[67..71].copy_from_slice(&2u32.to_be_bytes()); + section + } + fn build_mercator_section(scanning_mode: u8) -> Vec { let mut section = vec![0u8; 72]; section[..4].copy_from_slice(&72u32.to_be_bytes()); diff --git a/grib-core/src/lib.rs b/grib-core/src/lib.rs index ffde473..9a8114d 100644 --- a/grib-core/src/lib.rs +++ b/grib-core/src/lib.rs @@ -21,7 +21,7 @@ pub use data::{ pub use error::{Error, Result}; pub use grid::{ AlbersEqualAreaGrid, GridDefinition, LambertConformalGrid, LatLonGrid, MercatorGrid, - PolarStereographicGrid, ProjectedGridCore, + PolarStereographicGrid, ProjectedGridCore, RegularGaussianGrid, RotatedLatLonGrid, }; pub use metadata::{ForecastTimeUnit, Parameter, ParameterTableSource, ReferenceTime}; pub use parameter::{ diff --git a/grib-reader/src/grid.rs b/grib-reader/src/grid.rs index 245020c..60a95d3 100644 --- a/grib-reader/src/grid.rs +++ b/grib-reader/src/grid.rs @@ -1,4 +1,4 @@ pub use grib_core::grid::{ AlbersEqualAreaGrid, GridDefinition, LambertConformalGrid, LatLonGrid, MercatorGrid, - PolarStereographicGrid, ProjectedGridCore, + PolarStereographicGrid, ProjectedGridCore, RegularGaussianGrid, RotatedLatLonGrid, }; diff --git a/grib-reader/src/lib.rs b/grib-reader/src/lib.rs index a09ec2e..03bceb2 100644 --- a/grib-reader/src/lib.rs +++ b/grib-reader/src/lib.rs @@ -1,10 +1,9 @@ //! Pure-Rust GRIB file reader. //! -//! The current implementation supports the production-critical baseline for both -//! GRIB1 and GRIB2: regular latitude/longitude grids, GRIB2 Lambert conformal -//! and polar stereographic metadata and flat decode, simple packing, GRIB2 -//! complex packing with general group splitting, and optional image-backed -//! GRIB2 packing codecs. +//! The reader supports GRIB1 and GRIB2 regular geographic grids, GRIB2 rotated +//! latitude/longitude, regular Gaussian, and common projected grids, simple and +//! complex packing, multi-field bitmap reuse, bounded allocation, and optional +//! image-backed packing codecs. //! //! # Example //! @@ -51,6 +50,7 @@ pub use error::{Error, Result}; pub use grib1::{BinaryDataSection, GridDescription, ProductDefinition as Grib1ProductDefinition}; pub use grid::{ GridDefinition, LambertConformalGrid, LatLonGrid, PolarStereographicGrid, ProjectedGridCore, + RegularGaussianGrid, RotatedLatLonGrid, }; pub use metadata::{ForecastTimeUnit, Parameter, ParameterTableSource, ReferenceTime}; pub use parameter::{ @@ -504,16 +504,30 @@ impl<'a> Message<'a> { pub fn latitudes(&self) -> Result>> { self.metadata .grid - .as_lat_lon() - .map(|grid| grid.latitudes_with_limit(self.options.max_axis_points)) - .transpose() + .latitudes_with_limit(self.options.max_axis_points) } pub fn longitudes(&self) -> Result>> { self.metadata .grid - .as_lat_lon() - .map(|grid| grid.longitudes_with_limit(self.options.max_axis_points)) + .longitudes_with_limit(self.options.max_axis_points) + } + + /// Latitude axis in the rotated coordinate system for Template 3.1. + pub fn rotated_latitudes(&self) -> Result>> { + self.metadata + .grid + .as_rotated_lat_lon() + .map(|grid| grid.rotated_latitudes_with_limit(self.options.max_axis_points)) + .transpose() + } + + /// Longitude axis in the rotated coordinate system for Template 3.1. + pub fn rotated_longitudes(&self) -> Result>> { + self.metadata + .grid + .as_rotated_lat_lon() + .map(|grid| grid.rotated_longitudes_with_limit(self.options.max_axis_points)) .transpose() } @@ -1062,6 +1076,44 @@ mod tests { section } + fn build_rotated_grid() -> Vec { + let mut section = vec![0u8; 84]; + section[..4].copy_from_slice(&84u32.to_be_bytes()); + section[4] = 3; + section[6..10].copy_from_slice(&6u32.to_be_bytes()); + section[12..14].copy_from_slice(&1u16.to_be_bytes()); + section[30..34].copy_from_slice(&3u32.to_be_bytes()); + section[34..38].copy_from_slice(&2u32.to_be_bytes()); + section[46..50].copy_from_slice(&grib_i32_bytes(-10_000_000)); + section[50..54].copy_from_slice(&grib_i32_bytes(-20_000_000)); + section[55..59].copy_from_slice(&grib_i32_bytes(0)); + section[59..63].copy_from_slice(&grib_i32_bytes(20_000_000)); + section[63..67].copy_from_slice(&20_000_000u32.to_be_bytes()); + section[67..71].copy_from_slice(&10_000_000u32.to_be_bytes()); + section[71] = 0b0100_0000; + section[72..76].copy_from_slice(&grib_i32_bytes(-30_000_000)); + section[76..80].copy_from_slice(&10_000_000u32.to_be_bytes()); + section[80..84].copy_from_slice(&15.5f32.to_be_bytes()); + section + } + + fn build_regular_gaussian_grid() -> Vec { + let mut section = vec![0u8; 72]; + section[..4].copy_from_slice(&72u32.to_be_bytes()); + section[4] = 3; + section[6..10].copy_from_slice(&16u32.to_be_bytes()); + section[12..14].copy_from_slice(&40u16.to_be_bytes()); + section[30..34].copy_from_slice(&4u32.to_be_bytes()); + section[34..38].copy_from_slice(&4u32.to_be_bytes()); + section[46..50].copy_from_slice(&grib_i32_bytes(59_444_408)); + section[50..54].copy_from_slice(&grib_i32_bytes(0)); + section[55..59].copy_from_slice(&grib_i32_bytes(-59_444_408)); + section[59..63].copy_from_slice(&grib_i32_bytes(270_000_000)); + section[63..67].copy_from_slice(&90_000_000u32.to_be_bytes()); + section[67..71].copy_from_slice(&2u32.to_be_bytes()); + section + } + fn build_product(parameter_category: u8, parameter_number: u8) -> Vec { let mut section = vec![0u8; 34]; section[..4].copy_from_slice(&(34u32).to_be_bytes()); @@ -1260,6 +1312,63 @@ mod tests { ); } + #[test] + fn decodes_rotated_latlon_grid_and_exposes_rotated_axes() { + let message = assemble_grib2_message(&[ + build_identification(), + build_rotated_grid(), + build_product(0, 0), + build_simple_representation(6, 8), + build_data(&pack_u8_values(&[1, 2, 3, 4, 5, 6])), + ]); + let file = GribFile::from_bytes(message).unwrap(); + let field = file.message(0).unwrap(); + + assert_eq!(field.grid_shape(), (3, 2)); + assert_eq!(field.latitudes().unwrap(), None); + assert_eq!( + field.rotated_latitudes().unwrap().unwrap(), + vec![-10.0, 0.0] + ); + assert_eq!( + field.rotated_longitudes().unwrap().unwrap(), + vec![-20.0, 0.0, 20.0] + ); + assert_eq!( + field + .read_data_as_f64() + .unwrap() + .iter() + .copied() + .collect::>(), + vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0] + ); + } + + #[test] + fn decodes_regular_gaussian_grid_and_computes_axes() { + let values: Vec<_> = (0u8..16).collect(); + let message = assemble_grib2_message(&[ + build_identification(), + build_regular_gaussian_grid(), + build_product(0, 0), + build_simple_representation(values.len(), 8), + build_data(&pack_u8_values(&values)), + ]); + let file = GribFile::from_bytes(message).unwrap(); + let field = file.message(0).unwrap(); + + assert_eq!(field.grid_shape(), (4, 4)); + assert_eq!( + field.longitudes().unwrap().unwrap(), + vec![0.0, 90.0, 180.0, 270.0] + ); + let latitudes = field.latitudes().unwrap().unwrap(); + assert!((latitudes[0] - 59.444_408_289_166_77).abs() < 1e-12); + assert!((latitudes[3] + 59.444_408_289_166_77).abs() < 1e-12); + assert_eq!(field.read_data_as_f64().unwrap().shape(), &[4, 4]); + } + #[test] fn applies_bitmap_to_missing_values() { let message = assemble_grib2_message(&[ diff --git a/grib-reader/tests/common/fixtures.rs b/grib-reader/tests/common/fixtures.rs index 6c44075..84d2f44 100644 --- a/grib-reader/tests/common/fixtures.rs +++ b/grib-reader/tests/common/fixtures.rs @@ -93,6 +93,30 @@ pub fn build_grib2_polar_stereographic_alternating_message() -> Vec { build_grib2_polar_stereographic_message_with_scanning_mode(0b0001_0000, &[1, 2, 3, 6, 5, 4]) } +pub fn build_grib2_rotated_latlon_message() -> Vec { + let values = [1, 2, 3, 4, 5, 6]; + let sections = [ + build_identification_section(), + build_rotated_latlon_grid_section(), + build_product_section(), + build_simple_representation_section(values.len(), 8), + build_data_section(&values), + ]; + assemble_grib2_message(§ions) +} + +pub fn build_grib2_regular_gaussian_message() -> Vec { + let values: Vec<_> = (0u8..16).collect(); + let sections = [ + build_identification_section(), + build_regular_gaussian_grid_section(), + build_product_section(), + build_simple_representation_section(values.len(), 8), + build_data_section(&values), + ]; + assemble_grib2_message(§ions) +} + fn build_grib2_lambert_message_with_scanning_mode(scanning_mode: u8, values: &[u8]) -> Vec { let sections = [ build_identification_section(), @@ -180,6 +204,44 @@ fn build_polar_stereographic_grid_section(scanning_mode: u8) -> Vec { section } +fn build_rotated_latlon_grid_section() -> Vec { + let mut section = vec![0u8; 84]; + section[..4].copy_from_slice(&84u32.to_be_bytes()); + section[4] = 3; + section[6..10].copy_from_slice(&6u32.to_be_bytes()); + section[12..14].copy_from_slice(&1u16.to_be_bytes()); + section[30..34].copy_from_slice(&3u32.to_be_bytes()); + section[34..38].copy_from_slice(&2u32.to_be_bytes()); + section[46..50].copy_from_slice(&encode_wmo_i32(-10_000_000).unwrap()); + section[50..54].copy_from_slice(&encode_wmo_i32(-20_000_000).unwrap()); + section[55..59].copy_from_slice(&encode_wmo_i32(0).unwrap()); + section[59..63].copy_from_slice(&encode_wmo_i32(20_000_000).unwrap()); + section[63..67].copy_from_slice(&20_000_000u32.to_be_bytes()); + section[67..71].copy_from_slice(&10_000_000u32.to_be_bytes()); + section[71] = 0b0100_0000; + section[72..76].copy_from_slice(&encode_wmo_i32(-30_000_000).unwrap()); + section[76..80].copy_from_slice(&10_000_000u32.to_be_bytes()); + section[80..84].copy_from_slice(&15.5f32.to_be_bytes()); + section +} + +fn build_regular_gaussian_grid_section() -> Vec { + let mut section = vec![0u8; 72]; + section[..4].copy_from_slice(&72u32.to_be_bytes()); + section[4] = 3; + section[6..10].copy_from_slice(&16u32.to_be_bytes()); + section[12..14].copy_from_slice(&40u16.to_be_bytes()); + section[30..34].copy_from_slice(&4u32.to_be_bytes()); + section[34..38].copy_from_slice(&4u32.to_be_bytes()); + section[46..50].copy_from_slice(&encode_wmo_i32(59_444_408).unwrap()); + section[50..54].copy_from_slice(&encode_wmo_i32(0).unwrap()); + section[55..59].copy_from_slice(&encode_wmo_i32(-59_444_408).unwrap()); + section[59..63].copy_from_slice(&encode_wmo_i32(270_000_000).unwrap()); + section[63..67].copy_from_slice(&90_000_000u32.to_be_bytes()); + section[67..71].copy_from_slice(&2u32.to_be_bytes()); + section +} + fn build_product_section() -> Vec { let mut section = vec![0u8; 34]; section[..4].copy_from_slice(&34u32.to_be_bytes()); diff --git a/grib-reader/tests/common/mod.rs b/grib-reader/tests/common/mod.rs index ad759ec..715fa5f 100644 --- a/grib-reader/tests/common/mod.rs +++ b/grib-reader/tests/common/mod.rs @@ -9,7 +9,8 @@ pub use fixtures::{ build_grib2_complex_packing_message_with_missing, build_grib2_lambert_alternating_message, build_grib2_lambert_message, build_grib2_message, build_grib2_message_with_forecast, build_grib2_multifield_message, build_grib2_polar_stereographic_alternating_message, - build_grib2_polar_stereographic_message, build_grib2_spatial_differencing_message, + build_grib2_polar_stereographic_message, build_grib2_regular_gaussian_message, + build_grib2_rotated_latlon_message, build_grib2_spatial_differencing_message, build_truncated_grib2_message, }; diff --git a/grib-reader/tests/parity_eccodes.rs b/grib-reader/tests/parity_eccodes.rs index 6991aa0..d1325ba 100644 --- a/grib-reader/tests/parity_eccodes.rs +++ b/grib-reader/tests/parity_eccodes.rs @@ -6,6 +6,7 @@ use common::{ build_grib1_bitmap_message, build_grib1_message, build_grib2_complex_packing_message, build_grib2_complex_packing_message_with_missing, build_grib2_lambert_message, build_grib2_message, build_grib2_multifield_message, build_grib2_polar_stereographic_message, + build_grib2_regular_gaussian_message, build_grib2_rotated_latlon_message, build_grib2_spatial_differencing_message, collect_parity_samples, dump_reference, helper_path, write_fixture, }; @@ -57,6 +58,16 @@ fn generated_fixtures_match_eccodes_when_configured() { "polar-stereographic.grib2", &build_grib2_polar_stereographic_message(), ), + write_fixture( + dir.path(), + "rotated-latlon.grib2", + &build_grib2_rotated_latlon_message(), + ), + write_fixture( + dir.path(), + "regular-gaussian.grib2", + &build_grib2_regular_gaussian_message(), + ), ]; for path in fixtures { From b9cf6912078b0190071d4ab5c17f385b0e043420 Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 05:56:28 -0500 Subject: [PATCH 2/7] Preserve GRIB2 process metadata and signed forecast offsets --- CHANGELOG.md | 2 + README.md | 12 +- grib-core/src/metadata.rs | 34 +++++- grib-core/src/product.rs | 85 ++++++++++++-- .../fuzz_targets/fuzz_grib_writer_inputs.rs | 8 +- grib-reader/src/lib.rs | 11 +- grib-writer/src/lib.rs | 105 ++++++++++++++++-- grib-writer/tests/common/mod.rs | 29 ++++- 8 files changed, 252 insertions(+), 34 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 87f3eb6..4aef224 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,8 @@ offset - decode rotated latitude/longitude and regular Gaussian grids while rejecting quasi-regular geometry with a typed error +- preserve GRIB2 generating-process and data-cutoff metadata and decode signed + forecast offsets without unsigned wraparound - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index 0e12192..3a9c54d 100644 --- a/README.md +++ b/README.md @@ -132,14 +132,14 @@ let product = ProductDefinition { parameter_category: 0, parameter_number: 0, template: ProductDefinitionTemplate::AnalysisOrForecast(AnalysisOrForecastTemplate { - generating_process: 2, + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), forecast_time_unit: 1, forecast_time: 6, - first_surface: Some(FixedSurface { - surface_type: 103, - scale_factor: 0, - scaled_value: 850, - }), + first_surface: Some(FixedSurface::with_value(103, 0, 850)), second_surface: None, }), }; diff --git a/grib-core/src/metadata.rs b/grib-core/src/metadata.rs index f64bd1c..ca99f04 100644 --- a/grib-core/src/metadata.rs +++ b/grib-core/src/metadata.rs @@ -128,11 +128,11 @@ impl ReferenceTime { pub fn checked_add_forecast_time_unit( &self, unit: ForecastTimeUnit, - value: u32, + value: i64, ) -> Option { let seconds_per_unit = unit.seconds_per_unit()?; let base = self.seconds_since_epoch()?; - let delta = i64::from(value).checked_mul(seconds_per_unit)?; + let delta = value.checked_mul(seconds_per_unit)?; Self::from_seconds_since_epoch(base.checked_add(delta)?) } @@ -144,7 +144,7 @@ impl ReferenceTime { &self, edition: u8, unit: u8, - value: u32, + value: i64, ) -> Option { let unit = ForecastTimeUnit::from_edition_and_code(edition, unit)?; self.checked_add_forecast_time_unit(unit, value) @@ -154,7 +154,7 @@ impl ReferenceTime { /// /// Returns `None` for unsupported unit codes, calendar-dependent units, or /// invalid timestamps. - pub fn checked_add_forecast_time(&self, unit: u8, value: u32) -> Option { + pub fn checked_add_forecast_time(&self, unit: u8, value: i64) -> Option { let unit = ForecastTimeUnit::from_grib2_code(unit)?; self.checked_add_forecast_time_unit(unit, value) } @@ -371,6 +371,32 @@ mod tests { ); } + #[test] + fn subtracts_negative_forecast_offsets() { + let valid = ReferenceTime { + year: 2026, + month: 3, + day: 20, + hour: 6, + minute: 0, + second: 0, + } + .checked_add_forecast_time(1, -12) + .unwrap(); + + assert_eq!( + valid, + ReferenceTime { + year: 2026, + month: 3, + day: 19, + hour: 18, + minute: 0, + second: 0, + } + ); + } + #[test] fn adds_forecast_days_across_leap_day() { let valid = ReferenceTime { diff --git a/grib-core/src/product.rs b/grib-core/src/product.rs index 7721b1e..c5a64f6 100644 --- a/grib-core/src/product.rs +++ b/grib-core/src/product.rs @@ -135,9 +135,13 @@ pub enum ProductDefinitionTemplate { /// Product Definition Template 4.0: analysis or forecast at a horizontal level. #[derive(Debug, Clone, PartialEq)] pub struct AnalysisOrForecastTemplate { - pub generating_process: u8, + pub type_of_generating_process: u8, + pub background_generating_process_identifier: u8, + pub generating_process_identifier: u8, + pub hours_after_data_cutoff: Option, + pub minutes_after_data_cutoff: Option, pub forecast_time_unit: u8, - pub forecast_time: u32, + pub forecast_time: i32, pub first_surface: Option, pub second_surface: Option, } @@ -219,15 +223,23 @@ impl ProductDefinition { self.template.number() } - pub fn generating_process(&self) -> Option { - self.template.base().map(|base| base.generating_process) + pub fn type_of_generating_process(&self) -> Option { + self.template + .base() + .map(|base| base.type_of_generating_process) + } + + pub fn generating_process_identifier(&self) -> Option { + self.template + .base() + .map(|base| base.generating_process_identifier) } pub fn forecast_time_unit(&self) -> Option { self.template.base().map(|base| base.forecast_time_unit) } - pub fn forecast_time(&self) -> Option { + pub fn forecast_time(&self) -> Option { self.template.base().map(|base| base.forecast_time) } @@ -313,10 +325,23 @@ impl AnalysisOrForecastTemplate { fn parse(section_bytes: &[u8]) -> Result { require_len(section_bytes, Self::MINIMUM_LENGTH, "template 4.0")?; + let minutes_after_data_cutoff = (section_bytes[16] != 0xff).then_some(section_bytes[16]); + if minutes_after_data_cutoff.is_some_and(|minutes| minutes > 59) { + return Err(Error::InvalidSection { + section: 4, + reason: "minutes after data cutoff must be at most 59".into(), + }); + } + Ok(Self { - generating_process: section_bytes[11], + type_of_generating_process: section_bytes[11], + background_generating_process_identifier: section_bytes[12], + generating_process_identifier: section_bytes[13], + hours_after_data_cutoff: (section_bytes[14..16] != [0xff; 2]) + .then(|| u16::from_be_bytes(section_bytes[14..16].try_into().unwrap())), + minutes_after_data_cutoff, forecast_time_unit: section_bytes[17], - forecast_time: u32::from_be_bytes(section_bytes[18..22].try_into().unwrap()), + forecast_time: decode_wmo_i32(§ion_bytes[18..22]).unwrap(), first_surface: parse_surface(§ion_bytes[22..28]), second_surface: parse_surface(§ion_bytes[28..34]), }) @@ -511,7 +536,11 @@ mod tests { assert_eq!( product.template, ProductDefinitionTemplate::AnalysisOrForecast(AnalysisOrForecastTemplate { - generating_process: 2, + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), forecast_time_unit: 1, forecast_time: 6, first_surface: product.first_surface().cloned(), @@ -544,6 +573,46 @@ mod tests { } } + #[test] + fn parses_signed_forecast_time_and_process_metadata() { + let mut section = product_section_template_zero(); + section[12] = 7; + section[13] = 42; + section[14..16].copy_from_slice(&12u16.to_be_bytes()); + section[16] = 30; + section[18..22].copy_from_slice(&crate::binary::encode_wmo_i32(-6).unwrap()); + + let product = ProductDefinition::parse(§ion).unwrap(); + assert_eq!(product.type_of_generating_process(), Some(2)); + assert_eq!(product.generating_process_identifier(), Some(42)); + assert_eq!(product.forecast_time(), Some(-6)); + let ProductDefinitionTemplate::AnalysisOrForecast(template) = product.template else { + panic!("expected template 4.0"); + }; + assert_eq!(template.background_generating_process_identifier, 7); + assert_eq!(template.hours_after_data_cutoff, Some(12)); + assert_eq!(template.minutes_after_data_cutoff, Some(30)); + } + + #[test] + fn parses_missing_cutoff_metadata_and_rejects_invalid_minutes() { + let mut section = product_section_template_zero(); + section[14..16].copy_from_slice(&u16::MAX.to_be_bytes()); + section[16] = u8::MAX; + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::AnalysisOrForecast(template) = product.template else { + panic!("expected template 4.0"); + }; + assert_eq!(template.hours_after_data_cutoff, None); + assert_eq!(template.minutes_after_data_cutoff, None); + + section[16] = 60; + assert!(matches!( + ProductDefinition::parse(§ion), + Err(Error::InvalidSection { section: 4, .. }) + )); + } + #[test] fn parses_statistical_process_template() { let section = product_section_template_eight(); diff --git a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs index e4ad3eb..1df8cfd 100644 --- a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs +++ b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs @@ -217,9 +217,13 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { parameter_category, parameter_number, template: ProductDefinitionTemplate::AnalysisOrForecast(AnalysisOrForecastTemplate { - generating_process: 2, + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), forecast_time_unit: 1, - forecast_time: u32::from(input.u8()), + forecast_time: i32::from(input.u8()), first_surface: Some(FixedSurface::with_value(103, 0, 850)), second_surface: None, }), diff --git a/grib-reader/src/lib.rs b/grib-reader/src/lib.rs index 03bceb2..da96e49 100644 --- a/grib-reader/src/lib.rs +++ b/grib-reader/src/lib.rs @@ -235,7 +235,7 @@ pub struct MessageMetadata { pub grid: GridDefinition, pub data_representation: DataRepresentation, pub forecast_time_unit: Option, - pub forecast_time: Option, + pub forecast_time: Option, pub message_offset: u64, pub message_length: u64, pub field_index_in_message: usize, @@ -476,7 +476,10 @@ impl<'a> Message<'a> { ForecastTimeUnit::from_edition_and_code(self.metadata.edition, unit) } - pub fn forecast_time(&self) -> Option { + /// Forecast offset in the units returned by [`Message::forecast_time_unit`]. + /// + /// GRIB2 offsets are signed; GRIB1 offsets are widened to the same type. + pub fn forecast_time(&self) -> Option { self.metadata.forecast_time } @@ -890,7 +893,7 @@ fn index_grib1_message( grid, data_representation, forecast_time_unit: Some(sections.product.forecast_time_unit), - forecast_time: sections.product.forecast_time(), + forecast_time: sections.product.forecast_time().map(i64::from), message_offset: offset as u64, message_length: message_bytes.len() as u64, field_index_in_message: 0, @@ -994,7 +997,7 @@ fn index_grib2_message( grid, data_representation, forecast_time_unit: product.forecast_time_unit(), - forecast_time: product.forecast_time(), + forecast_time: product.forecast_time().map(i64::from), message_offset: offset as u64, message_length: message_bytes.len() as u64, field_index_in_message, diff --git a/grib-writer/src/lib.rs b/grib-writer/src/lib.rs index aaf6275..0ddc634 100644 --- a/grib-writer/src/lib.rs +++ b/grib-writer/src/lib.rs @@ -2034,23 +2034,49 @@ fn write_product_template_prefix( section_length: u32, template: &AnalysisOrForecastTemplate, ) -> Result<()> { + validate_product_template_prefix(template)?; + write_u32_be(out, section_length)?; write_u8_be(out, 4)?; write_u16_be(out, 0)?; write_u16_be(out, template_number)?; write_u8_be(out, product.parameter_category)?; write_u8_be(out, product.parameter_number)?; - write_u8_be(out, template.generating_process)?; - write_u8_be(out, 0)?; - write_u8_be(out, 0)?; - write_u16_be(out, 0)?; - write_u8_be(out, 0)?; + write_u8_be(out, template.type_of_generating_process)?; + write_u8_be(out, template.background_generating_process_identifier)?; + write_u8_be(out, template.generating_process_identifier)?; + write_u16_be(out, template.hours_after_data_cutoff.unwrap_or(u16::MAX))?; + write_u8_be(out, template.minutes_after_data_cutoff.unwrap_or(u8::MAX))?; write_u8_be(out, template.forecast_time_unit)?; - write_u32_be(out, template.forecast_time)?; + out.extend_from_slice(&encode_wmo_i32(template.forecast_time).ok_or_else(|| { + Error::ValueOutOfRange("forecast time does not fit GRIB signed i32".into()) + })?); write_surface(out, template.first_surface.as_ref())?; write_surface(out, template.second_surface.as_ref()) } +fn validate_product_template_prefix(template: &AnalysisOrForecastTemplate) -> Result<()> { + if template.hours_after_data_cutoff == Some(u16::MAX) { + return Err(Error::ValueOutOfRange( + "hours after data cutoff must be at most 65534".into(), + )); + } + if template + .minutes_after_data_cutoff + .is_some_and(|minutes| minutes > 59) + { + return Err(Error::ValueOutOfRange( + "minutes after data cutoff must be at most 59".into(), + )); + } + if encode_wmo_i32(template.forecast_time).is_none() { + return Err(Error::ValueOutOfRange( + "forecast time does not fit GRIB signed i32".into(), + )); + } + Ok(()) +} + fn write_ensemble_product_extra( out: &mut Vec, template: &grib_core::IndividualEnsembleForecastTemplate, @@ -2397,13 +2423,19 @@ fn validate_supported_grib1_grid(grid: &GridDefinition) -> Result<()> { fn validate_supported_product(product: &ProductDefinition) -> Result<()> { match &product.template { - ProductDefinitionTemplate::AnalysisOrForecast(_) => Ok(()), - ProductDefinitionTemplate::IndividualEnsembleForecast(_) => Ok(()), + ProductDefinitionTemplate::AnalysisOrForecast(template) => { + validate_product_template_prefix(template) + } + ProductDefinitionTemplate::IndividualEnsembleForecast(template) => { + validate_product_template_prefix(&template.base) + } ProductDefinitionTemplate::StatisticalProcess(template) => { + validate_product_template_prefix(&template.base)?; checked_time_range_count(template.time_ranges.len())?; validate_reference_time(template.end_of_overall_time_interval) } ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { + validate_product_template_prefix(&template.ensemble.base)?; checked_time_range_count(template.time_ranges.len())?; validate_reference_time(template.end_of_overall_time_interval) } @@ -2641,7 +2673,11 @@ mod tests { fn analysis_or_forecast_template() -> AnalysisOrForecastTemplate { AnalysisOrForecastTemplate { - generating_process: 2, + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), forecast_time_unit: 1, forecast_time: 6, first_surface: Some(FixedSurface::with_value(103, 0, 850)), @@ -2909,6 +2945,57 @@ mod tests { assert_eq!(message.read_flat_data_as_f64().unwrap(), values); } + #[test] + fn writes_signed_forecast_offset_and_process_metadata() { + let mut template = analysis_or_forecast_template(); + template.background_generating_process_identifier = 7; + template.generating_process_identifier = 42; + template.hours_after_data_cutoff = Some(12); + template.minutes_after_data_cutoff = Some(30); + template.forecast_time = -6; + let field = Grib2FieldBuilder::new() + .identification(identification()) + .grid(grid()) + .product(ProductDefinition { + parameter_category: 0, + parameter_number: 0, + template: ProductDefinitionTemplate::AnalysisOrForecast(template), + }) + .packing(PackingStrategy::SimpleAuto { decimal_scale: 0 }) + .values(&[1.0, 2.0, 3.0, 4.0]) + .build() + .unwrap(); + + let file = GribFile::from_bytes(write_message([field])).unwrap(); + let message = file.message(0).unwrap(); + assert_eq!(message.forecast_time(), Some(-6)); + let product = message.product_definition().unwrap(); + assert_eq!(product.generating_process_identifier(), Some(42)); + let ProductDefinitionTemplate::AnalysisOrForecast(template) = &product.template else { + panic!("expected template 4.0"); + }; + assert_eq!(template.background_generating_process_identifier, 7); + assert_eq!(template.hours_after_data_cutoff, Some(12)); + assert_eq!(template.minutes_after_data_cutoff, Some(30)); + } + + #[test] + fn rejects_unrepresentable_cutoff_metadata() { + let mut template = analysis_or_forecast_template(); + template.hours_after_data_cutoff = Some(u16::MAX); + let mut product = product(0, 0); + product.template = ProductDefinitionTemplate::AnalysisOrForecast(template); + let err = Grib2FieldBuilder::new() + .identification(identification()) + .grid(grid()) + .product(product) + .packing(PackingStrategy::SimpleAuto { decimal_scale: 0 }) + .values(&[1.0, 2.0, 3.0, 4.0]) + .build() + .unwrap_err(); + assert!(matches!(err, grib_core::Error::ValueOutOfRange(_))); + } + #[test] fn writes_individual_ensemble_product_template_readable_by_reader() { let values = [1.0, 2.0, 3.0, 4.0]; diff --git a/grib-writer/tests/common/mod.rs b/grib-writer/tests/common/mod.rs index a0f0075..94bbefe 100644 --- a/grib-writer/tests/common/mod.rs +++ b/grib-writer/tests/common/mod.rs @@ -244,6 +244,25 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { .values(&spatial_second_values) .build() .unwrap(); + let mut signed_forecast_product = product(0, 0); + let ProductDefinitionTemplate::AnalysisOrForecast(template) = + &mut signed_forecast_product.template + else { + unreachable!("product helper always returns template 4.0"); + }; + template.background_generating_process_identifier = 7; + template.generating_process_identifier = 42; + template.hours_after_data_cutoff = Some(12); + template.minutes_after_data_cutoff = Some(30); + template.forecast_time = -6; + let signed_forecast = Grib2FieldBuilder::new() + .identification(identification()) + .grid(latlon_grid(2, 2, 0)) + .product(signed_forecast_product) + .packing(PackingStrategy::SimpleAuto { decimal_scale: 0 }) + .values(&[1.0, 2.0, 3.0, 4.0]) + .build() + .unwrap(); vec![ ( @@ -255,6 +274,10 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { write_grib2_message([simple_grib2_field(&[5.0, f64::NAN, 7.0, 8.0], 0, 0)]), ), ("writer-decimal.grib2", write_grib2_message([decimal])), + ( + "writer-signed-forecast.grib2", + write_grib2_message([signed_forecast]), + ), ("writer-complex.grib2", write_grib2_message([complex])), ( "writer-complex-spatial-first.grib2", @@ -393,7 +416,11 @@ pub fn product(parameter_category: u8, parameter_number: u8) -> ProductDefinitio parameter_category, parameter_number, template: ProductDefinitionTemplate::AnalysisOrForecast(AnalysisOrForecastTemplate { - generating_process: 2, + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), forecast_time_unit: 1, forecast_time: 6, first_surface: Some(FixedSurface::with_value(103, 0, 850)), From 8ea3ae7b2db988f3bb541878c4713750189f0152 Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 07:47:43 -0500 Subject: [PATCH 3/7] Decode derived probability and percentile product templates --- CHANGELOG.md | 2 + README.md | 12 +- grib-core/src/lib.rs | 7 +- grib-core/src/product.rs | 329 ++++++++++++++++++++++++++++++++++++- grib-reader/src/lib.rs | 7 +- grib-reader/src/product.rs | 7 +- 6 files changed, 350 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4aef224..89531f5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,8 @@ quasi-regular geometry with a typed error - preserve GRIB2 generating-process and data-cutoff metadata and decode signed forecast offsets without unsigned wraparound +- decode derived-ensemble, probability, and percentile forecasts with typed + probability thresholds for product templates 4.2, 4.5, and 4.6 - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index 3a9c54d..f0f0ef3 100644 --- a/README.md +++ b/README.md @@ -176,7 +176,9 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Typed metadata access for reference time, parameter identity, product metadata, grid geometry, and lat/lon coordinates - Explicit north-up decode methods in addition to reader order that preserves the encoded i/j scan directions -- Reader and writer GRIB2 product definition templates 4.0, 4.1, 4.8, and 4.11 +- Reader GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, + and 4.11 +- Writer GRIB2 product definition templates 4.0, 4.1, 4.8, and 4.11 - Forecast valid-time helpers for supported fixed-width GRIB1/GRIB2 time units - `GribFile::builder()` for strict or tolerant scanning and allocation limits - Bitmap application with missing values surfaced as `NaN` @@ -198,10 +200,14 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Remaining GRIB2 grid templates beyond 3.0, 3.1, 3.10, 3.20, 3.30, 3.31, and 3.40 -- Remaining GRIB2 product definition templates beyond 4.0, 4.1, 4.8, and 4.11 +- Remaining reader GRIB2 product definition templates beyond 4.0, 4.1, 4.2, + 4.5, 4.6, 4.8, and 4.11 +- Remaining writer GRIB2 product definition templates beyond 4.0, 4.1, 4.8, + and 4.11 - Writer GRIB2 row-by-row complex packing -Unsupported cases fail explicitly with typed errors. +Unsupported decode and encode operations fail explicitly with typed errors; +unknown GRIB2 product templates remain indexable with their raw bytes preserved. Calendar-dependent forecast units such as months and years are exposed through raw metadata but currently return `None` from `valid_time()`. diff --git a/grib-core/src/lib.rs b/grib-core/src/lib.rs index 9a8114d..9a6ce7d 100644 --- a/grib-core/src/lib.rs +++ b/grib-core/src/lib.rs @@ -29,7 +29,8 @@ pub use parameter::{ LOCAL_PARAMETER_TABLE_CSV_HEADER, }; pub use product::{ - AnalysisOrForecastTemplate, EnsembleStatisticalProcessTemplate, FixedSurface, Identification, - IndividualEnsembleForecastTemplate, ProductDefinition, ProductDefinitionTemplate, ScaledValue, - StatisticalProcessTemplate, StatisticalTimeRange, + AnalysisOrForecastTemplate, DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, + FixedSurface, Identification, IndividualEnsembleForecastTemplate, PercentileForecastTemplate, + ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ScaledValue, StatisticalProcessTemplate, StatisticalTimeRange, }; diff --git a/grib-core/src/product.rs b/grib-core/src/product.rs index c5a64f6..5325171 100644 --- a/grib-core/src/product.rs +++ b/grib-core/src/product.rs @@ -121,6 +121,9 @@ pub struct ProductDefinition { pub enum ProductDefinitionTemplate { AnalysisOrForecast(AnalysisOrForecastTemplate), IndividualEnsembleForecast(IndividualEnsembleForecastTemplate), + DerivedForecast(DerivedForecastTemplate), + ProbabilityForecast(ProbabilityForecastTemplate), + PercentileForecast(PercentileForecastTemplate), StatisticalProcess(StatisticalProcessTemplate), EnsembleStatisticalProcess(EnsembleStatisticalProcessTemplate), /// A well-framed Section 4 whose template is not interpreted by this @@ -155,6 +158,136 @@ pub struct IndividualEnsembleForecastTemplate { pub number_of_forecasts_in_ensemble: u8, } +/// Product Definition Template 4.2: forecast derived from all ensemble members. +#[derive(Debug, Clone, PartialEq)] +pub struct DerivedForecastTemplate { + pub base: AnalysisOrForecastTemplate, + pub derived_forecast_type: u8, + pub number_of_forecasts_in_ensemble: u8, +} + +/// A signed decimal threshold used by probability product templates. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct ProbabilityLimit { + pub scale_factor: i16, + pub scaled_value: i32, +} + +impl ProbabilityLimit { + pub fn value_f64(self) -> f64 { + f64::from(self.scaled_value) * 10.0_f64.powi(-i32::from(self.scale_factor)) + } +} + +/// The event whose forecast probability is encoded by templates 4.5 and 4.9. +/// +/// Known WMO event types carry exactly the limit values used by their +/// definition. [`ProbabilityType::Other`] preserves reserved, locally defined, +/// and noncanonical input without imposing semantics on its limits. The writer +/// accepts `Other` only for codes without a standard typed variant. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[non_exhaustive] +pub enum ProbabilityType { + BelowLowerLimit(ProbabilityLimit), + AboveUpperLimit(ProbabilityLimit), + BetweenLimits { + lower: ProbabilityLimit, + upper: ProbabilityLimit, + }, + AboveLowerLimit(ProbabilityLimit), + BelowUpperLimit(ProbabilityLimit), + EqualLowerLimit(ProbabilityLimit), + AboveNormal, + NearNormal, + BelowNormal, + CategoricalBoolean, + Quantile, + Missing, + Other { + code: u8, + lower: Option, + upper: Option, + }, +} + +impl ProbabilityType { + pub const fn code(self) -> u8 { + match self { + Self::BelowLowerLimit(_) => 0, + Self::AboveUpperLimit(_) => 1, + Self::BetweenLimits { .. } => 2, + Self::AboveLowerLimit(_) => 3, + Self::BelowUpperLimit(_) => 4, + Self::EqualLowerLimit(_) => 5, + Self::AboveNormal => 6, + Self::NearNormal => 7, + Self::BelowNormal => 8, + Self::CategoricalBoolean => 9, + Self::Quantile => 10, + Self::Missing => 255, + Self::Other { code, .. } => code, + } + } + + pub const fn lower_limit(self) -> Option { + match self { + Self::BelowLowerLimit(limit) + | Self::AboveLowerLimit(limit) + | Self::EqualLowerLimit(limit) => Some(limit), + Self::BetweenLimits { lower, .. } => Some(lower), + Self::Other { lower, .. } => lower, + _ => None, + } + } + + pub const fn upper_limit(self) -> Option { + match self { + Self::AboveUpperLimit(limit) | Self::BelowUpperLimit(limit) => Some(limit), + Self::BetweenLimits { upper, .. } => Some(upper), + Self::Other { upper, .. } => upper, + _ => None, + } + } + + fn from_code_and_limits( + code: u8, + lower: Option, + upper: Option, + ) -> Self { + match (code, lower, upper) { + (0, Some(limit), None) => Self::BelowLowerLimit(limit), + (1, None, Some(limit)) => Self::AboveUpperLimit(limit), + (2, Some(lower), Some(upper)) => Self::BetweenLimits { lower, upper }, + (3, Some(limit), None) => Self::AboveLowerLimit(limit), + (4, None, Some(limit)) => Self::BelowUpperLimit(limit), + (5, Some(limit), None) => Self::EqualLowerLimit(limit), + (6, None, None) => Self::AboveNormal, + (7, None, None) => Self::NearNormal, + (8, None, None) => Self::BelowNormal, + (9, None, None) => Self::CategoricalBoolean, + (10, None, None) => Self::Quantile, + (255, None, None) => Self::Missing, + (code, lower, upper) => Self::Other { code, lower, upper }, + } + } +} + +/// Product Definition Template 4.5: probability forecast at a point in time. +#[derive(Debug, Clone, PartialEq)] +pub struct ProbabilityForecastTemplate { + pub base: AnalysisOrForecastTemplate, + pub forecast_probability_number: u8, + pub total_number_of_forecast_probabilities: u8, + pub probability: ProbabilityType, +} + +/// Product Definition Template 4.6: percentile forecast at a point in time. +#[derive(Debug, Clone, PartialEq)] +pub struct PercentileForecastTemplate { + pub base: AnalysisOrForecastTemplate, + pub percentile_value: u8, +} + /// Product Definition Template 4.8: statistically processed field over a time interval. #[derive(Debug, Clone, PartialEq)] pub struct StatisticalProcessTemplate { @@ -269,6 +402,15 @@ impl ProductDefinitionTemplate { 1 => Ok(Self::IndividualEnsembleForecast( IndividualEnsembleForecastTemplate::parse(section_bytes)?, )), + 2 => Ok(Self::DerivedForecast(DerivedForecastTemplate::parse( + section_bytes, + )?)), + 5 => Ok(Self::ProbabilityForecast( + ProbabilityForecastTemplate::parse(section_bytes)?, + )), + 6 => Ok(Self::PercentileForecast(PercentileForecastTemplate::parse( + section_bytes, + )?)), 8 => Ok(Self::StatisticalProcess(StatisticalProcessTemplate::parse( section_bytes, )?)), @@ -291,6 +433,9 @@ impl ProductDefinitionTemplate { match self { Self::AnalysisOrForecast(_) => 0, Self::IndividualEnsembleForecast(_) => 1, + Self::DerivedForecast(_) => 2, + Self::ProbabilityForecast(_) => 5, + Self::PercentileForecast(_) => 6, Self::StatisticalProcess(_) => 8, Self::EnsembleStatisticalProcess(_) => 11, Self::Unsupported { number, .. } => *number, @@ -301,6 +446,9 @@ impl ProductDefinitionTemplate { Some(match self { Self::AnalysisOrForecast(template) => template, Self::IndividualEnsembleForecast(template) => &template.base, + Self::DerivedForecast(template) => &template.base, + Self::ProbabilityForecast(template) => &template.base, + Self::PercentileForecast(template) => &template.base, Self::StatisticalProcess(template) => &template.base, Self::EnsembleStatisticalProcess(template) => &template.ensemble.base, Self::Unsupported { .. } => return None, @@ -313,7 +461,11 @@ impl ProductDefinitionTemplate { Self::EnsembleStatisticalProcess(template) => { Some(template.end_of_overall_time_interval) } - Self::AnalysisOrForecast(_) | Self::IndividualEnsembleForecast(_) => None, + Self::AnalysisOrForecast(_) + | Self::IndividualEnsembleForecast(_) + | Self::DerivedForecast(_) + | Self::ProbabilityForecast(_) + | Self::PercentileForecast(_) => None, Self::Unsupported { .. } => None, } } @@ -363,6 +515,57 @@ impl IndividualEnsembleForecastTemplate { } } +impl DerivedForecastTemplate { + const MINIMUM_LENGTH: usize = 36; + + fn parse(section_bytes: &[u8]) -> Result { + require_len(section_bytes, Self::MINIMUM_LENGTH, "template 4.2")?; + + Ok(Self { + base: AnalysisOrForecastTemplate::parse(section_bytes)?, + derived_forecast_type: section_bytes[34], + number_of_forecasts_in_ensemble: section_bytes[35], + }) + } +} + +impl ProbabilityForecastTemplate { + const MINIMUM_LENGTH: usize = 47; + + fn parse(section_bytes: &[u8]) -> Result { + require_len(section_bytes, Self::MINIMUM_LENGTH, "template 4.5")?; + + let lower = parse_probability_limit(§ion_bytes[37..42]); + let upper = parse_probability_limit(§ion_bytes[42..47]); + Ok(Self { + base: AnalysisOrForecastTemplate::parse(section_bytes)?, + forecast_probability_number: section_bytes[34], + total_number_of_forecast_probabilities: section_bytes[35], + probability: ProbabilityType::from_code_and_limits(section_bytes[36], lower, upper), + }) + } +} + +impl PercentileForecastTemplate { + const MINIMUM_LENGTH: usize = 35; + + fn parse(section_bytes: &[u8]) -> Result { + require_len(section_bytes, Self::MINIMUM_LENGTH, "template 4.6")?; + + let percentile_value = section_bytes[34]; + if percentile_value > 100 { + return Err(Error::InvalidSection { + section: 4, + reason: format!("template 4.6 percentile {percentile_value} exceeds 100"), + }); + } + Ok(Self { + base: AnalysisOrForecastTemplate::parse(section_bytes)?, + percentile_value, + }) + } +} + impl StatisticalProcessTemplate { const TIME_RANGE_OFFSET: usize = 46; @@ -466,6 +669,17 @@ fn parse_statistical_time_ranges( .collect() } +fn parse_probability_limit(bytes: &[u8]) -> Option { + if bytes[0] == 0xff || bytes[1..5] == [0xff; 4] { + return None; + } + + Some(ProbabilityLimit { + scale_factor: decode_wmo_i8(bytes[0]), + scaled_value: decode_wmo_i32(&bytes[1..5])?, + }) +} + fn parse_surface(section_bytes: &[u8]) -> Option { let surface_type = section_bytes[0]; if surface_type == 255 { @@ -490,7 +704,8 @@ fn parse_surface(section_bytes: &[u8]) -> Option { #[cfg(test)] mod tests { use super::{ - AnalysisOrForecastTemplate, Identification, ProductDefinition, ProductDefinitionTemplate, + AnalysisOrForecastTemplate, Identification, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, }; use crate::error::Error; use crate::metadata::ReferenceTime; @@ -573,6 +788,97 @@ mod tests { } } + #[test] + fn parses_derived_forecast_template() { + let mut section = product_section_template_zero(); + section.resize(36, 0); + set_product_template(&mut section, 2); + section[34] = 4; + section[35] = 50; + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::DerivedForecast(template) = product.template else { + panic!("expected template 4.2"); + }; + assert_eq!(template.derived_forecast_type, 4); + assert_eq!(template.number_of_forecasts_in_ensemble, 50); + assert_eq!(template.base.forecast_time, 6); + } + + #[test] + fn parses_probability_forecast_with_typed_thresholds() { + let mut section = product_section_template_zero(); + section.resize(47, 0xff); + set_product_template(&mut section, 5); + section[34] = 2; + section[35] = 10; + section[36] = 2; + section[37] = 1; + section[38..42].copy_from_slice(&crate::binary::encode_wmo_i32(-125).unwrap()); + section[42] = 1; + section[43..47].copy_from_slice(&crate::binary::encode_wmo_i32(250).unwrap()); + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::ProbabilityForecast(template) = product.template else { + panic!("expected template 4.5"); + }; + assert_eq!(template.forecast_probability_number, 2); + assert_eq!(template.total_number_of_forecast_probabilities, 10); + let ProbabilityType::BetweenLimits { lower, upper } = template.probability else { + panic!("expected a between-limits probability"); + }; + assert_eq!(lower.value_f64(), -12.5); + assert_eq!(upper.value_f64(), 25.0); + } + + #[test] + fn preserves_noncanonical_probability_limit_combinations() { + let mut section = product_section_template_zero(); + section.resize(47, 0xff); + set_product_template(&mut section, 5); + section[36] = 0; + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::ProbabilityForecast(template) = product.template else { + panic!("expected template 4.5"); + }; + assert_eq!( + template.probability, + ProbabilityType::Other { + code: 0, + lower: None, + upper: None, + } + ); + + section[36] = 7; + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::ProbabilityForecast(template) = product.template else { + panic!("expected template 4.5"); + }; + assert_eq!(template.probability, ProbabilityType::NearNormal); + } + + #[test] + fn parses_percentile_forecast_and_rejects_values_above_one_hundred() { + let mut section = product_section_template_zero(); + section.resize(35, 0); + set_product_template(&mut section, 6); + section[34] = 90; + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::PercentileForecast(template) = product.template else { + panic!("expected template 4.6"); + }; + assert_eq!(template.percentile_value, 90); + + section[34] = 101; + assert!(matches!( + ProductDefinition::parse(§ion), + Err(Error::InvalidSection { section: 4, .. }) + )); + } + #[test] fn parses_signed_forecast_time_and_process_metadata() { let mut section = product_section_template_zero(); @@ -733,6 +1039,19 @@ mod tests { assert!(matches!(err, Error::InvalidSection { section: 4, .. })); } + #[test] + fn rejects_truncated_instantaneous_product_templates() { + for (template, length) in [(2, 35), (5, 46), (6, 34)] { + let mut section = product_section_template_zero(); + section.resize(length, 0); + set_product_template(&mut section, template); + assert!(matches!( + ProductDefinition::parse(§ion), + Err(Error::InvalidSection { section: 4, .. }) + )); + } + } + fn product_section_template_zero() -> Vec { let mut section = vec![0u8; 34]; section[..4].copy_from_slice(&(34u32).to_be_bytes()); @@ -770,6 +1089,12 @@ mod tests { section } + fn set_product_template(section: &mut [u8], template: u16) { + let length = u32::try_from(section.len()).unwrap(); + section[..4].copy_from_slice(&length.to_be_bytes()); + section[7..9].copy_from_slice(&template.to_be_bytes()); + } + fn valid_identification_section() -> Vec { let mut section = vec![0u8; 21]; section[..4].copy_from_slice(&(21u32).to_be_bytes()); diff --git a/grib-reader/src/lib.rs b/grib-reader/src/lib.rs index da96e49..dcda18d 100644 --- a/grib-reader/src/lib.rs +++ b/grib-reader/src/lib.rs @@ -58,9 +58,10 @@ pub use parameter::{ LOCAL_PARAMETER_TABLE_CSV_HEADER, }; pub use product::{ - AnalysisOrForecastTemplate, EnsembleStatisticalProcessTemplate, FixedSurface, Identification, - IndividualEnsembleForecastTemplate, ProductDefinition, ProductDefinitionTemplate, ScaledValue, - StatisticalProcessTemplate, StatisticalTimeRange, + AnalysisOrForecastTemplate, DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, + FixedSurface, Identification, IndividualEnsembleForecastTemplate, PercentileForecastTemplate, + ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ScaledValue, StatisticalProcessTemplate, StatisticalTimeRange, }; use std::io::Read; diff --git a/grib-reader/src/product.rs b/grib-reader/src/product.rs index de85cd8..f8ef688 100644 --- a/grib-reader/src/product.rs +++ b/grib-reader/src/product.rs @@ -1,5 +1,6 @@ pub use grib_core::product::{ - AnalysisOrForecastTemplate, EnsembleStatisticalProcessTemplate, FixedSurface, Identification, - IndividualEnsembleForecastTemplate, ProductDefinition, ProductDefinitionTemplate, ScaledValue, - StatisticalProcessTemplate, StatisticalTimeRange, + AnalysisOrForecastTemplate, DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, + FixedSurface, Identification, IndividualEnsembleForecastTemplate, PercentileForecastTemplate, + ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ScaledValue, StatisticalProcessTemplate, StatisticalTimeRange, }; From 620489c01dcc5ef297339a9033a9a81f45c26057 Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 10:24:33 -0500 Subject: [PATCH 4/7] Write product templates 4.2 4.5 and 4.6 and verify ecCodes metadata --- CHANGELOG.md | 2 + README.md | 7 +- .../fuzz_targets/fuzz_grib_writer_inputs.rs | 68 ++++- grib-writer/src/lib.rs | 266 +++++++++++++++++- grib-writer/tests/common/mod.rs | 172 ++++++++++- tools/eccodes-reference.c | 125 ++++++++ 6 files changed, 605 insertions(+), 35 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 89531f5..bf39d7b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,8 @@ forecast offsets without unsigned wraparound - decode derived-ensemble, probability, and percentile forecasts with typed probability thresholds for product templates 4.2, 4.5, and 4.6 +- write product templates 4.2, 4.5, and 4.6 and compare their Section 4 + metadata with ecCodes in the parity suite - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index f0f0ef3..30c28c4 100644 --- a/README.md +++ b/README.md @@ -178,7 +178,8 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; the encoded i/j scan directions - Reader GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, and 4.11 -- Writer GRIB2 product definition templates 4.0, 4.1, 4.8, and 4.11 +- Writer GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, + and 4.11 - Forecast valid-time helpers for supported fixed-width GRIB1/GRIB2 time units - `GribFile::builder()` for strict or tolerant scanning and allocation limits - Bitmap application with missing values surfaced as `NaN` @@ -202,8 +203,8 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; and 3.40 - Remaining reader GRIB2 product definition templates beyond 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, and 4.11 -- Remaining writer GRIB2 product definition templates beyond 4.0, 4.1, 4.8, - and 4.11 +- Remaining writer GRIB2 product definition templates beyond 4.0, 4.1, 4.2, + 4.5, 4.6, 4.8, and 4.11 - Writer GRIB2 row-by-row complex packing Unsupported decode and encode operations fail explicitly with typed errors; diff --git a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs index 1df8cfd..1c692af 100644 --- a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs +++ b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs @@ -2,8 +2,9 @@ use grib_core::metadata::ReferenceTime; use grib_core::{ - AnalysisOrForecastTemplate, FixedSurface, GridDefinition, Identification, LatLonGrid, - ProductDefinition, ProductDefinitionTemplate, + AnalysisOrForecastTemplate, DerivedForecastTemplate, FixedSurface, GridDefinition, + Identification, LatLonGrid, PercentileForecastTemplate, ProbabilityForecastTemplate, + ProbabilityLimit, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, }; use grib_reader::GribFile; use grib_writer::{ @@ -213,20 +214,61 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { 2 => (0, 4), _ => (input.u8(), input.u8()), }; + let base = AnalysisOrForecastTemplate { + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), + forecast_time_unit: 1, + forecast_time: i32::from(input.u8()), + first_surface: Some(FixedSurface::with_value(103, 0, 850)), + second_surface: None, + }; + let template = match input.u8() % 4 { + 0 => ProductDefinitionTemplate::AnalysisOrForecast(base), + 1 => ProductDefinitionTemplate::DerivedForecast(DerivedForecastTemplate { + base, + derived_forecast_type: input.u8(), + number_of_forecasts_in_ensemble: input.u8(), + }), + 2 => ProductDefinitionTemplate::ProbabilityForecast(ProbabilityForecastTemplate { + base, + forecast_probability_number: input.u8(), + total_number_of_forecast_probabilities: input.u8(), + probability: probability_type(input), + }), + _ => ProductDefinitionTemplate::PercentileForecast(PercentileForecastTemplate { + base, + percentile_value: input.u8() % 101, + }), + }; ProductDefinition { parameter_category, parameter_number, - template: ProductDefinitionTemplate::AnalysisOrForecast(AnalysisOrForecastTemplate { - type_of_generating_process: 2, - background_generating_process_identifier: 0, - generating_process_identifier: 0, - hours_after_data_cutoff: Some(0), - minutes_after_data_cutoff: Some(0), - forecast_time_unit: 1, - forecast_time: i32::from(input.u8()), - first_surface: Some(FixedSurface::with_value(103, 0, 850)), - second_surface: None, - }), + template, + } +} + +fn probability_type(input: &mut Input<'_>) -> ProbabilityType { + let limit = |input: &mut Input<'_>| ProbabilityLimit { + scale_factor: i16::from(input.u8() % 15) - 7, + scaled_value: i32::from(input.i16()), + }; + match input.u8() % 6 { + 0 => ProbabilityType::BelowLowerLimit(limit(input)), + 1 => ProbabilityType::AboveUpperLimit(limit(input)), + 2 => ProbabilityType::BetweenLimits { + lower: limit(input), + upper: limit(input), + }, + 3 => ProbabilityType::NearNormal, + 4 => ProbabilityType::Quantile, + _ => ProbabilityType::Other { + code: 192 + input.u8() % 63, + lower: input.bool().then(|| limit(input)), + upper: input.bool().then(|| limit(input)), + }, } } diff --git a/grib-writer/src/lib.rs b/grib-writer/src/lib.rs index 0ddc634..ef162f3 100644 --- a/grib-writer/src/lib.rs +++ b/grib-writer/src/lib.rs @@ -13,8 +13,9 @@ use grib_core::{ AlbersEqualAreaGrid, AnalysisOrForecastTemplate, ComplexPackingParams, DataRepresentation, FixedSurface, GridDefinition, Identification, ImagePackingParams, Jpeg2000PackingParams, LambertConformalGrid, LatLonGrid, MercatorGrid, PngPackingParams, PolarStereographicGrid, - ProductDefinition, ProductDefinitionTemplate, ProjectedGridCore, ReferenceTime, - SimplePackingParams, SpatialDifferencingParams, StatisticalTimeRange, + ProbabilityLimit, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, + ProjectedGridCore, ReferenceTime, SimplePackingParams, SpatialDifferencingParams, + StatisticalTimeRange, }; pub use grib_core::grib1::ProductDefinition as Grib1ProductDefinition; @@ -1995,6 +1996,21 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu write_product_template_prefix(out, product, 1, 37, &template.base)?; write_ensemble_product_extra(out, template) } + ProductDefinitionTemplate::DerivedForecast(template) => { + write_product_template_prefix(out, product, 2, 36, &template.base)?; + write_u8_be(out, template.derived_forecast_type)?; + write_u8_be(out, template.number_of_forecasts_in_ensemble) + } + ProductDefinitionTemplate::ProbabilityForecast(template) => { + validate_probability_type(template.probability)?; + write_product_template_prefix(out, product, 5, 47, &template.base)?; + write_probability_product_extra(out, template) + } + ProductDefinitionTemplate::PercentileForecast(template) => { + validate_percentile(template.percentile_value)?; + write_product_template_prefix(out, product, 6, 35, &template.base)?; + write_u8_be(out, template.percentile_value) + } ProductDefinitionTemplate::StatisticalProcess(template) => { let range_count = checked_time_range_count(template.time_ranges.len())?; let section_length = statistical_product_section_len(46, range_count)?; @@ -2086,6 +2102,74 @@ fn write_ensemble_product_extra( write_u8_be(out, template.number_of_forecasts_in_ensemble) } +fn write_probability_product_extra( + out: &mut Vec, + template: &grib_core::ProbabilityForecastTemplate, +) -> Result<()> { + write_u8_be(out, template.forecast_probability_number)?; + write_u8_be(out, template.total_number_of_forecast_probabilities)?; + write_u8_be(out, template.probability.code())?; + write_probability_limit(out, template.probability.lower_limit())?; + write_probability_limit(out, template.probability.upper_limit()) +} + +fn write_probability_limit(out: &mut Vec, limit: Option) -> Result<()> { + let Some(limit) = limit else { + out.extend_from_slice(&[0xff; 5]); + return Ok(()); + }; + + write_u8_be( + out, + encode_wmo_i8(limit.scale_factor).ok_or_else(|| { + Error::ValueOutOfRange( + "probability-limit scale factor does not fit GRIB signed i8".into(), + ) + })?, + )?; + out.extend_from_slice(&encode_wmo_i32(limit.scaled_value).ok_or_else(|| { + Error::ValueOutOfRange("probability-limit scaled value does not fit GRIB signed i32".into()) + })?); + Ok(()) +} + +fn validate_probability_type(probability: ProbabilityType) -> Result<()> { + if let ProbabilityType::Other { code, .. } = probability { + if matches!(code, 0..=10 | 255) { + return Err(Error::ValueOutOfRange(format!( + "WMO probability type {code} must use its typed probability variant" + ))); + } + } + + for limit in [probability.lower_limit(), probability.upper_limit()] + .into_iter() + .flatten() + { + if encode_wmo_i8(limit.scale_factor).is_none() { + return Err(Error::ValueOutOfRange( + "probability-limit scale factor does not fit GRIB signed i8".into(), + )); + } + if encode_wmo_i32(limit.scaled_value).is_none() { + return Err(Error::ValueOutOfRange( + "probability-limit scaled value does not fit GRIB signed i32".into(), + )); + } + } + Ok(()) +} + +fn validate_percentile(percentile: u8) -> Result<()> { + if percentile <= 100 { + Ok(()) + } else { + Err(Error::ValueOutOfRange(format!( + "percentile value {percentile} exceeds 100" + ))) + } +} + fn write_reference_time(out: &mut Vec, reference_time: ReferenceTime) -> Result<()> { validate_reference_time(reference_time)?; @@ -2429,6 +2513,17 @@ fn validate_supported_product(product: &ProductDefinition) -> Result<()> { ProductDefinitionTemplate::IndividualEnsembleForecast(template) => { validate_product_template_prefix(&template.base) } + ProductDefinitionTemplate::DerivedForecast(template) => { + validate_product_template_prefix(&template.base) + } + ProductDefinitionTemplate::ProbabilityForecast(template) => { + validate_product_template_prefix(&template.base)?; + validate_probability_type(template.probability) + } + ProductDefinitionTemplate::PercentileForecast(template) => { + validate_product_template_prefix(&template.base)?; + validate_percentile(template.percentile_value) + } ProductDefinitionTemplate::StatisticalProcess(template) => { validate_product_template_prefix(&template.base)?; checked_time_range_count(template.time_ranges.len())?; @@ -2458,10 +2553,12 @@ mod tests { use grib_core::metadata::ReferenceTime; use grib_core::{ AlbersEqualAreaGrid, AnalysisOrForecastTemplate, DataRepresentation, - EnsembleStatisticalProcessTemplate, FixedSurface, GridDefinition, Identification, - IndividualEnsembleForecastTemplate, LambertConformalGrid, LatLonGrid, MercatorGrid, - PolarStereographicGrid, ProductDefinition, ProductDefinitionTemplate, ProjectedGridCore, - StatisticalProcessTemplate, StatisticalTimeRange, + DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, FixedSurface, GridDefinition, + Identification, IndividualEnsembleForecastTemplate, LambertConformalGrid, LatLonGrid, + MercatorGrid, PercentileForecastTemplate, PolarStereographicGrid, + ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ProjectedGridCore, StatisticalProcessTemplate, + StatisticalTimeRange, }; use grib_reader::sections::scan_sections; use grib_reader::{GribFile, PredefinedBitmap}; @@ -2756,6 +2853,21 @@ mod tests { .unwrap() } + fn field_with_product_template(template: ProductDefinitionTemplate) -> super::Grib2Field { + Grib2FieldBuilder::new() + .identification(identification()) + .grid(grid()) + .product(ProductDefinition { + parameter_category: 0, + parameter_number: 0, + template, + }) + .packing(PackingStrategy::SimpleAuto { decimal_scale: 0 }) + .values(&[1.0, 2.0, 3.0, 4.0]) + .build() + .unwrap() + } + fn grib1_simple_field(values: &[f64]) -> super::Grib1Field { Grib1FieldBuilder::new() .product(grib1_product()) @@ -3043,6 +3155,148 @@ mod tests { assert_eq!(message.read_flat_data_as_f64().unwrap(), values); } + #[test] + fn writes_derived_forecast_product_template_readable_by_reader() { + let field = field_with_product_template(ProductDefinitionTemplate::DerivedForecast( + DerivedForecastTemplate { + base: analysis_or_forecast_template(), + derived_forecast_type: 4, + number_of_forecasts_in_ensemble: 50, + }, + )); + + let bytes = write_message([field]); + let product_section = scan_sections(&bytes) + .unwrap() + .into_iter() + .find(|section| section.number == 4) + .unwrap(); + assert_eq!(product_section.length, 36); + + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::DerivedForecast(template) = &product.template else { + panic!("expected template 4.2"); + }; + assert_eq!(template.derived_forecast_type, 4); + assert_eq!(template.number_of_forecasts_in_ensemble, 50); + assert_eq!( + message.read_flat_data_as_f64().unwrap(), + [1.0, 2.0, 3.0, 4.0] + ); + } + + #[test] + fn writes_probability_forecast_product_template_readable_by_reader() { + let probability = ProbabilityType::BetweenLimits { + lower: ProbabilityLimit { + scale_factor: 1, + scaled_value: -125, + }, + upper: ProbabilityLimit { + scale_factor: 1, + scaled_value: 250, + }, + }; + let field = field_with_product_template(ProductDefinitionTemplate::ProbabilityForecast( + ProbabilityForecastTemplate { + base: analysis_or_forecast_template(), + forecast_probability_number: 2, + total_number_of_forecast_probabilities: 10, + probability, + }, + )); + + let bytes = write_message([field]); + let product_section = scan_sections(&bytes) + .unwrap() + .into_iter() + .find(|section| section.number == 4) + .unwrap(); + assert_eq!(product_section.length, 47); + + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::ProbabilityForecast(template) = &product.template else { + panic!("expected template 4.5"); + }; + assert_eq!(template.forecast_probability_number, 2); + assert_eq!(template.total_number_of_forecast_probabilities, 10); + assert_eq!(template.probability, probability); + assert_eq!( + message.read_flat_data_as_f64().unwrap(), + [1.0, 2.0, 3.0, 4.0] + ); + } + + #[test] + fn writes_percentile_forecast_product_template_readable_by_reader() { + let field = field_with_product_template(ProductDefinitionTemplate::PercentileForecast( + PercentileForecastTemplate { + base: analysis_or_forecast_template(), + percentile_value: 90, + }, + )); + + let bytes = write_message([field]); + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::PercentileForecast(template) = &product.template else { + panic!("expected template 4.6"); + }; + assert_eq!(template.percentile_value, 90); + assert_eq!( + message.read_flat_data_as_f64().unwrap(), + [1.0, 2.0, 3.0, 4.0] + ); + } + + #[test] + fn rejects_noncanonical_probability_types_and_invalid_percentiles() { + for template in [ + ProductDefinitionTemplate::ProbabilityForecast(ProbabilityForecastTemplate { + base: analysis_or_forecast_template(), + forecast_probability_number: 1, + total_number_of_forecast_probabilities: 10, + probability: ProbabilityType::Other { + code: 0, + lower: None, + upper: None, + }, + }), + ProductDefinitionTemplate::ProbabilityForecast(ProbabilityForecastTemplate { + base: analysis_or_forecast_template(), + forecast_probability_number: 1, + total_number_of_forecast_probabilities: 10, + probability: ProbabilityType::BelowLowerLimit(ProbabilityLimit { + scale_factor: 128, + scaled_value: 10, + }), + }), + ProductDefinitionTemplate::PercentileForecast(PercentileForecastTemplate { + base: analysis_or_forecast_template(), + percentile_value: 101, + }), + ] { + let err = Grib2FieldBuilder::new() + .identification(identification()) + .grid(grid()) + .product(ProductDefinition { + parameter_category: 0, + parameter_number: 0, + template, + }) + .packing(PackingStrategy::SimpleAuto { decimal_scale: 0 }) + .values(&[1.0, 2.0, 3.0, 4.0]) + .build() + .unwrap_err(); + assert!(matches!(err, grib_core::Error::ValueOutOfRange(_))); + } + } + #[test] fn writes_statistical_product_template_readable_by_reader() { let values = [1.0, 2.0, 3.0, 4.0]; diff --git a/grib-writer/tests/common/mod.rs b/grib-writer/tests/common/mod.rs index 94bbefe..14bd490 100644 --- a/grib-writer/tests/common/mod.rs +++ b/grib-writer/tests/common/mod.rs @@ -5,8 +5,9 @@ use std::process::Command; use grib_core::metadata::ReferenceTime; use grib_core::{ - AnalysisOrForecastTemplate, FixedSurface, GridDefinition, Identification, LatLonGrid, - ProductDefinition, ProductDefinitionTemplate, + AnalysisOrForecastTemplate, DerivedForecastTemplate, FixedSurface, GridDefinition, + Identification, LatLonGrid, PercentileForecastTemplate, ProbabilityForecastTemplate, + ProbabilityLimit, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, }; use grib_reader::GribFile; use grib_writer::{ @@ -27,6 +28,17 @@ pub struct ReferenceMessage { pub reference_time: ReferenceTimeDump, pub ni: usize, pub nj: usize, + pub product_definition_template_number: Option, + pub derived_forecast: Option, + pub number_of_forecasts_in_ensemble: Option, + pub forecast_probability_number: Option, + pub total_number_of_forecast_probabilities: Option, + pub probability_type: Option, + pub scale_factor_of_lower_limit: Option, + pub scaled_value_of_lower_limit: Option, + pub scale_factor_of_upper_limit: Option, + pub scaled_value_of_upper_limit: Option, + pub percentile_value: Option, pub values: Vec>, } @@ -150,6 +162,7 @@ pub fn assert_matches_reference(helper: &Path, path: &Path, bytes: &[u8]) { path.display(), index ); + assert_product_metadata(&message, expected, path, index); assert_eq!( actual.len(), expected.values.len(), @@ -190,6 +203,91 @@ pub fn assert_matches_reference(helper: &Path, path: &Path, bytes: &[u8]) { } } +fn assert_product_metadata( + message: &grib_reader::Message<'_>, + expected: &ReferenceMessage, + path: &Path, + field_index: usize, +) { + if message.edition() != 2 { + return; + } + + let product = message.product_definition().unwrap_or_else(|| { + panic!( + "missing product metadata for {} field {}", + path.display(), + field_index + ) + }); + assert_eq!( + expected.product_definition_template_number, + Some(i64::from(product.template_number())), + "product template mismatch for {} field {}", + path.display(), + field_index + ); + + match &product.template { + ProductDefinitionTemplate::DerivedForecast(template) => { + assert_eq!( + expected.derived_forecast, + Some(i64::from(template.derived_forecast_type)) + ); + assert_eq!( + expected.number_of_forecasts_in_ensemble, + Some(i64::from(template.number_of_forecasts_in_ensemble)) + ); + } + ProductDefinitionTemplate::ProbabilityForecast(template) => { + assert_eq!( + expected.forecast_probability_number, + Some(i64::from(template.forecast_probability_number)) + ); + assert_eq!( + expected.total_number_of_forecast_probabilities, + Some(i64::from(template.total_number_of_forecast_probabilities)) + ); + assert_eq!( + expected.probability_type, + Some(i64::from(template.probability.code())) + ); + assert_probability_limit_metadata( + template.probability.lower_limit(), + expected.scale_factor_of_lower_limit, + expected.scaled_value_of_lower_limit, + ); + assert_probability_limit_metadata( + template.probability.upper_limit(), + expected.scale_factor_of_upper_limit, + expected.scaled_value_of_upper_limit, + ); + } + ProductDefinitionTemplate::PercentileForecast(template) => { + assert_eq!( + expected.percentile_value, + Some(i64::from(template.percentile_value)) + ); + } + _ => {} + } +} + +fn assert_probability_limit_metadata( + actual: Option, + expected_scale_factor: Option, + expected_scaled_value: Option, +) { + assert_eq!( + actual.map(|limit| i64::from(limit.scale_factor)), + expected_scale_factor + ); + assert_eq!( + actual.map(|limit| i64::from(limit.scaled_value)), + expected_scaled_value + ); +} + pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { let decimal = Grib2FieldBuilder::new() .identification(identification()) @@ -263,6 +361,44 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { .values(&[1.0, 2.0, 3.0, 4.0]) .build() .unwrap(); + let product_field = |template| { + Grib2FieldBuilder::new() + .identification(identification()) + .grid(latlon_grid(2, 2, 0)) + .product(ProductDefinition { + parameter_category: 0, + parameter_number: 0, + template, + }) + .packing(PackingStrategy::SimpleAuto { decimal_scale: 0 }) + .values(&[1.0, 2.0, 3.0, 4.0]) + .build() + .unwrap() + }; + let derived = product_field(ProductDefinitionTemplate::DerivedForecast( + DerivedForecastTemplate { + base: analysis_or_forecast_template(), + derived_forecast_type: 4, + number_of_forecasts_in_ensemble: 50, + }, + )); + let probability = product_field(ProductDefinitionTemplate::ProbabilityForecast( + ProbabilityForecastTemplate { + base: analysis_or_forecast_template(), + forecast_probability_number: 1, + total_number_of_forecast_probabilities: 10, + probability: ProbabilityType::BelowLowerLimit(ProbabilityLimit { + scale_factor: 1, + scaled_value: 2732, + }), + }, + )); + let percentile = product_field(ProductDefinitionTemplate::PercentileForecast( + PercentileForecastTemplate { + base: analysis_or_forecast_template(), + percentile_value: 90, + }, + )); vec![ ( @@ -278,6 +414,12 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { "writer-signed-forecast.grib2", write_grib2_message([signed_forecast]), ), + ("writer-derived.grib2", write_grib2_message([derived])), + ( + "writer-probability.grib2", + write_grib2_message([probability]), + ), + ("writer-percentile.grib2", write_grib2_message([percentile])), ("writer-complex.grib2", write_grib2_message([complex])), ( "writer-complex-spatial-first.grib2", @@ -415,17 +557,21 @@ pub fn product(parameter_category: u8, parameter_number: u8) -> ProductDefinitio ProductDefinition { parameter_category, parameter_number, - template: ProductDefinitionTemplate::AnalysisOrForecast(AnalysisOrForecastTemplate { - type_of_generating_process: 2, - background_generating_process_identifier: 0, - generating_process_identifier: 0, - hours_after_data_cutoff: Some(0), - minutes_after_data_cutoff: Some(0), - forecast_time_unit: 1, - forecast_time: 6, - first_surface: Some(FixedSurface::with_value(103, 0, 850)), - second_surface: None, - }), + template: ProductDefinitionTemplate::AnalysisOrForecast(analysis_or_forecast_template()), + } +} + +pub fn analysis_or_forecast_template() -> AnalysisOrForecastTemplate { + AnalysisOrForecastTemplate { + type_of_generating_process: 2, + background_generating_process_identifier: 0, + generating_process_identifier: 0, + hours_after_data_cutoff: Some(0), + minutes_after_data_cutoff: Some(0), + forecast_time_unit: 1, + forecast_time: 6, + first_surface: Some(FixedSurface::with_value(103, 0, 850)), + second_surface: None, } } diff --git a/tools/eccodes-reference.c b/tools/eccodes-reference.c index fdbf3d0..18b3fbf 100644 --- a/tools/eccodes-reference.c +++ b/tools/eccodes-reference.c @@ -55,6 +55,27 @@ static long get_long_or_default( return value; } +static int get_optional_long( + codes_handle *handle, + const char *key, + long *value, + const char *path +) { + int err = CODES_SUCCESS; + int missing = codes_is_missing(handle, key, &err); + if (err == CODES_NOT_FOUND) { + return 0; + } + if (err != CODES_SUCCESS) { + die_codes(err, key, path); + } + if (missing) { + return 0; + } + *value = get_long(handle, key, path); + return 1; +} + static long get_grid_dimension( codes_handle *handle, const char *primary_key, @@ -249,6 +270,14 @@ static void print_json_string(const char *value) { fputc('"', stdout); } +static void print_optional_long(int present, long value) { + if (present) { + fprintf(stdout, "%ld", value); + } else { + fputs("null", stdout); + } +} + static decode_totals decode_file(const char *path, int emit_json) { FILE *fp = fopen(path, "rb"); if (fp == NULL) { @@ -280,6 +309,71 @@ static decode_totals decode_file(const char *path, int emit_json) { long second = get_long(handle, "second", path); long ni = get_grid_dimension(handle, "Ni", "Nx", path); long nj = get_grid_dimension(handle, "Nj", "Ny", path); + long product_definition_template_number = 0; + long derived_forecast = 0; + long number_of_forecasts_in_ensemble = 0; + long forecast_probability_number = 0; + long total_number_of_forecast_probabilities = 0; + long probability_type = 0; + long scale_factor_of_lower_limit = 0; + long scaled_value_of_lower_limit = 0; + long scale_factor_of_upper_limit = 0; + long scaled_value_of_upper_limit = 0; + long percentile_value = 0; + int has_product_definition_template_number = get_optional_long( + handle, + "productDefinitionTemplateNumber", + &product_definition_template_number, + path + ); + int has_derived_forecast = + get_optional_long(handle, "derivedForecast", &derived_forecast, path); + int has_number_of_forecasts_in_ensemble = get_optional_long( + handle, + "numberOfForecastsInEnsemble", + &number_of_forecasts_in_ensemble, + path + ); + int has_forecast_probability_number = get_optional_long( + handle, + "forecastProbabilityNumber", + &forecast_probability_number, + path + ); + int has_total_number_of_forecast_probabilities = get_optional_long( + handle, + "totalNumberOfForecastProbabilities", + &total_number_of_forecast_probabilities, + path + ); + int has_probability_type = + get_optional_long(handle, "probabilityType", &probability_type, path); + int has_scale_factor_of_lower_limit = get_optional_long( + handle, + "scaleFactorOfLowerLimit", + &scale_factor_of_lower_limit, + path + ); + int has_scaled_value_of_lower_limit = get_optional_long( + handle, + "scaledValueOfLowerLimit", + &scaled_value_of_lower_limit, + path + ); + int has_scale_factor_of_upper_limit = get_optional_long( + handle, + "scaleFactorOfUpperLimit", + &scale_factor_of_upper_limit, + path + ); + int has_scaled_value_of_upper_limit = get_optional_long( + handle, + "scaledValueOfUpperLimit", + &scaled_value_of_upper_limit, + path + ); + int has_percentile_value = + get_optional_long(handle, "percentileValue", &percentile_value, path); char name[256]; get_string(handle, "name", name, sizeof(name), path); @@ -318,6 +412,37 @@ static decode_totals decode_file(const char *path, int emit_json) { fprintf(stdout, "%ld", ni); fputs(",\"nj\":", stdout); fprintf(stdout, "%ld", nj); + fputs(",\"product_definition_template_number\":", stdout); + print_optional_long( + has_product_definition_template_number, + product_definition_template_number + ); + fputs(",\"derived_forecast\":", stdout); + print_optional_long(has_derived_forecast, derived_forecast); + fputs(",\"number_of_forecasts_in_ensemble\":", stdout); + print_optional_long( + has_number_of_forecasts_in_ensemble, + number_of_forecasts_in_ensemble + ); + fputs(",\"forecast_probability_number\":", stdout); + print_optional_long(has_forecast_probability_number, forecast_probability_number); + fputs(",\"total_number_of_forecast_probabilities\":", stdout); + print_optional_long( + has_total_number_of_forecast_probabilities, + total_number_of_forecast_probabilities + ); + fputs(",\"probability_type\":", stdout); + print_optional_long(has_probability_type, probability_type); + fputs(",\"scale_factor_of_lower_limit\":", stdout); + print_optional_long(has_scale_factor_of_lower_limit, scale_factor_of_lower_limit); + fputs(",\"scaled_value_of_lower_limit\":", stdout); + print_optional_long(has_scaled_value_of_lower_limit, scaled_value_of_lower_limit); + fputs(",\"scale_factor_of_upper_limit\":", stdout); + print_optional_long(has_scale_factor_of_upper_limit, scale_factor_of_upper_limit); + fputs(",\"scaled_value_of_upper_limit\":", stdout); + print_optional_long(has_scaled_value_of_upper_limit, scaled_value_of_upper_limit); + fputs(",\"percentile_value\":", stdout); + print_optional_long(has_percentile_value, percentile_value); fputs(",\"values\":[", stdout); for (size_t i = 0; i < value_len; ++i) { if (i > 0) { From cf379be3791671d4b9e68165d011fc33d6565c82 Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 11:42:53 -0500 Subject: [PATCH 5/7] Decode product templates 4.9 4.10 and 4.12 through one interval model --- CHANGELOG.md | 2 + README.md | 4 +- grib-core/src/lib.rs | 11 +- grib-core/src/product.rs | 294 +++++++++++++++++++++++++------ grib-reader/src/lib.rs | 11 +- grib-reader/src/product.rs | 11 +- grib-reader/tests/integration.rs | 11 +- grib-writer/src/lib.rs | 90 ++++++---- 8 files changed, 320 insertions(+), 114 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bf39d7b..58d87a0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,8 @@ probability thresholds for product templates 4.2, 4.5, and 4.6 - write product templates 4.2, 4.5, and 4.6 and compare their Section 4 metadata with ecCodes in the parity suite +- decode interval probability, percentile, and derived-ensemble templates 4.9, + 4.10, and 4.12 through one shared statistical-interval model - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index 30c28c4..9b1a554 100644 --- a/README.md +++ b/README.md @@ -177,7 +177,7 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Explicit north-up decode methods in addition to reader order that preserves the encoded i/j scan directions - Reader GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, - and 4.11 + 4.9, 4.10, 4.11, and 4.12 - Writer GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, and 4.11 - Forecast valid-time helpers for supported fixed-width GRIB1/GRIB2 time units @@ -202,7 +202,7 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Remaining GRIB2 grid templates beyond 3.0, 3.1, 3.10, 3.20, 3.30, 3.31, and 3.40 - Remaining reader GRIB2 product definition templates beyond 4.0, 4.1, 4.2, - 4.5, 4.6, 4.8, and 4.11 + 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, and 4.12 - Remaining writer GRIB2 product definition templates beyond 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, and 4.11 - Writer GRIB2 row-by-row complex packing diff --git a/grib-core/src/lib.rs b/grib-core/src/lib.rs index 9a6ce7d..0ce23c5 100644 --- a/grib-core/src/lib.rs +++ b/grib-core/src/lib.rs @@ -29,8 +29,11 @@ pub use parameter::{ LOCAL_PARAMETER_TABLE_CSV_HEADER, }; pub use product::{ - AnalysisOrForecastTemplate, DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, - FixedSurface, Identification, IndividualEnsembleForecastTemplate, PercentileForecastTemplate, - ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ScaledValue, StatisticalProcessTemplate, StatisticalTimeRange, + AnalysisOrForecastTemplate, DerivedForecastTemplate, DerivedStatisticalProcessTemplate, + EnsembleStatisticalProcessTemplate, FixedSurface, Identification, + IndividualEnsembleForecastTemplate, PercentileForecastTemplate, + PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, + ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ScaledValue, StatisticalInterval, StatisticalProcessTemplate, + StatisticalTimeRange, }; diff --git a/grib-core/src/product.rs b/grib-core/src/product.rs index 5325171..750763e 100644 --- a/grib-core/src/product.rs +++ b/grib-core/src/product.rs @@ -125,7 +125,10 @@ pub enum ProductDefinitionTemplate { ProbabilityForecast(ProbabilityForecastTemplate), PercentileForecast(PercentileForecastTemplate), StatisticalProcess(StatisticalProcessTemplate), + ProbabilityStatisticalProcess(ProbabilityStatisticalProcessTemplate), + PercentileStatisticalProcess(PercentileStatisticalProcessTemplate), EnsembleStatisticalProcess(EnsembleStatisticalProcessTemplate), + DerivedStatisticalProcess(DerivedStatisticalProcessTemplate), /// A well-framed Section 4 whose template is not interpreted by this /// version of the library. `raw` contains the template-specific bytes /// following the common parameter category and number. @@ -288,22 +291,47 @@ pub struct PercentileForecastTemplate { pub percentile_value: u8, } -/// Product Definition Template 4.8: statistically processed field over a time interval. +/// Time-interval metadata shared by statistically processed product templates. #[derive(Debug, Clone, PartialEq)] -pub struct StatisticalProcessTemplate { - pub base: AnalysisOrForecastTemplate, +pub struct StatisticalInterval { pub end_of_overall_time_interval: ReferenceTime, pub number_of_missing_in_statistical_process: u32, pub time_ranges: Vec, } +/// Product Definition Template 4.8: statistically processed field over a time interval. +#[derive(Debug, Clone, PartialEq)] +pub struct StatisticalProcessTemplate { + pub base: AnalysisOrForecastTemplate, + pub interval: StatisticalInterval, +} + +/// Product Definition Template 4.9: probability forecast over a time interval. +#[derive(Debug, Clone, PartialEq)] +pub struct ProbabilityStatisticalProcessTemplate { + pub probability: ProbabilityForecastTemplate, + pub interval: StatisticalInterval, +} + +/// Product Definition Template 4.10: percentile forecast over a time interval. +#[derive(Debug, Clone, PartialEq)] +pub struct PercentileStatisticalProcessTemplate { + pub percentile: PercentileForecastTemplate, + pub interval: StatisticalInterval, +} + /// Product Definition Template 4.11: individual ensemble forecast over a time interval. #[derive(Debug, Clone, PartialEq)] pub struct EnsembleStatisticalProcessTemplate { pub ensemble: IndividualEnsembleForecastTemplate, - pub end_of_overall_time_interval: ReferenceTime, - pub number_of_missing_in_statistical_process: u32, - pub time_ranges: Vec, + pub interval: StatisticalInterval, +} + +/// Product Definition Template 4.12: derived ensemble forecast over a time interval. +#[derive(Debug, Clone, PartialEq)] +pub struct DerivedStatisticalProcessTemplate { + pub derived: DerivedForecastTemplate, + pub interval: StatisticalInterval, } /// Statistical processing descriptor from GRIB2 Product Definition templates @@ -414,9 +442,18 @@ impl ProductDefinitionTemplate { 8 => Ok(Self::StatisticalProcess(StatisticalProcessTemplate::parse( section_bytes, )?)), + 9 => Ok(Self::ProbabilityStatisticalProcess( + ProbabilityStatisticalProcessTemplate::parse(section_bytes)?, + )), + 10 => Ok(Self::PercentileStatisticalProcess( + PercentileStatisticalProcessTemplate::parse(section_bytes)?, + )), 11 => Ok(Self::EnsembleStatisticalProcess( EnsembleStatisticalProcessTemplate::parse(section_bytes)?, )), + 12 => Ok(Self::DerivedStatisticalProcess( + DerivedStatisticalProcessTemplate::parse(section_bytes)?, + )), number => { let raw_len = section_bytes.len() - 11; let mut raw = Vec::new(); @@ -437,7 +474,10 @@ impl ProductDefinitionTemplate { Self::ProbabilityForecast(_) => 5, Self::PercentileForecast(_) => 6, Self::StatisticalProcess(_) => 8, + Self::ProbabilityStatisticalProcess(_) => 9, + Self::PercentileStatisticalProcess(_) => 10, Self::EnsembleStatisticalProcess(_) => 11, + Self::DerivedStatisticalProcess(_) => 12, Self::Unsupported { number, .. } => *number, } } @@ -450,16 +490,30 @@ impl ProductDefinitionTemplate { Self::ProbabilityForecast(template) => &template.base, Self::PercentileForecast(template) => &template.base, Self::StatisticalProcess(template) => &template.base, + Self::ProbabilityStatisticalProcess(template) => &template.probability.base, + Self::PercentileStatisticalProcess(template) => &template.percentile.base, Self::EnsembleStatisticalProcess(template) => &template.ensemble.base, + Self::DerivedStatisticalProcess(template) => &template.derived.base, Self::Unsupported { .. } => return None, }) } fn end_of_overall_time_interval(&self) -> Option { match self { - Self::StatisticalProcess(template) => Some(template.end_of_overall_time_interval), + Self::StatisticalProcess(template) => { + Some(template.interval.end_of_overall_time_interval) + } + Self::ProbabilityStatisticalProcess(template) => { + Some(template.interval.end_of_overall_time_interval) + } + Self::PercentileStatisticalProcess(template) => { + Some(template.interval.end_of_overall_time_interval) + } Self::EnsembleStatisticalProcess(template) => { - Some(template.end_of_overall_time_interval) + Some(template.interval.end_of_overall_time_interval) + } + Self::DerivedStatisticalProcess(template) => { + Some(template.interval.end_of_overall_time_interval) } Self::AnalysisOrForecast(_) | Self::IndividualEnsembleForecast(_) @@ -567,47 +621,78 @@ impl PercentileForecastTemplate { } impl StatisticalProcessTemplate { - const TIME_RANGE_OFFSET: usize = 46; + fn parse(section_bytes: &[u8]) -> Result { + Ok(Self { + base: AnalysisOrForecastTemplate::parse(section_bytes)?, + interval: StatisticalInterval::parse(section_bytes, 34, "template 4.8")?, + }) + } +} +impl ProbabilityStatisticalProcessTemplate { fn parse(section_bytes: &[u8]) -> Result { - require_len(section_bytes, Self::TIME_RANGE_OFFSET, "template 4.8")?; - let time_range_count = section_bytes[41] as usize; - let min_len = required_time_range_template_len(Self::TIME_RANGE_OFFSET, time_range_count)?; - require_len(section_bytes, min_len, "template 4.8")?; + Ok(Self { + probability: ProbabilityForecastTemplate::parse(section_bytes)?, + interval: StatisticalInterval::parse(section_bytes, 47, "template 4.9")?, + }) + } +} +impl PercentileStatisticalProcessTemplate { + fn parse(section_bytes: &[u8]) -> Result { Ok(Self { - base: AnalysisOrForecastTemplate::parse(section_bytes)?, - end_of_overall_time_interval: parse_reference_time(§ion_bytes[34..41], 4)?, - number_of_missing_in_statistical_process: u32::from_be_bytes( - section_bytes[42..46].try_into().unwrap(), - ), - time_ranges: parse_statistical_time_ranges( - §ion_bytes[Self::TIME_RANGE_OFFSET..min_len], - time_range_count, - ), + percentile: PercentileForecastTemplate::parse(section_bytes)?, + interval: StatisticalInterval::parse(section_bytes, 35, "template 4.10")?, }) } } impl EnsembleStatisticalProcessTemplate { - const TIME_RANGE_OFFSET: usize = 49; + fn parse(section_bytes: &[u8]) -> Result { + Ok(Self { + ensemble: IndividualEnsembleForecastTemplate::parse(section_bytes)?, + interval: StatisticalInterval::parse(section_bytes, 37, "template 4.11")?, + }) + } +} +impl DerivedStatisticalProcessTemplate { fn parse(section_bytes: &[u8]) -> Result { - require_len(section_bytes, Self::TIME_RANGE_OFFSET, "template 4.11")?; - let time_range_count = section_bytes[44] as usize; - let min_len = required_time_range_template_len(Self::TIME_RANGE_OFFSET, time_range_count)?; - require_len(section_bytes, min_len, "template 4.11")?; + Ok(Self { + derived: DerivedForecastTemplate::parse(section_bytes)?, + interval: StatisticalInterval::parse(section_bytes, 36, "template 4.12")?, + }) + } +} + +impl StatisticalInterval { + fn parse(section_bytes: &[u8], end_time_offset: usize, context: &str) -> Result { + let time_range_offset = + end_time_offset + .checked_add(12) + .ok_or_else(|| Error::InvalidSection { + section: 4, + reason: "statistical interval offset overflow".into(), + })?; + require_len(section_bytes, time_range_offset, context)?; + let time_range_count = section_bytes[end_time_offset + 7] as usize; + let min_len = required_time_range_template_len(time_range_offset, time_range_count)?; + require_len(section_bytes, min_len, context)?; Ok(Self { - ensemble: IndividualEnsembleForecastTemplate::parse(section_bytes)?, - end_of_overall_time_interval: parse_reference_time(§ion_bytes[37..44], 4)?, + end_of_overall_time_interval: parse_reference_time( + §ion_bytes[end_time_offset..end_time_offset + 7], + 4, + )?, number_of_missing_in_statistical_process: u32::from_be_bytes( - section_bytes[45..49].try_into().unwrap(), + section_bytes[end_time_offset + 8..end_time_offset + 12] + .try_into() + .unwrap(), ), time_ranges: parse_statistical_time_ranges( - §ion_bytes[Self::TIME_RANGE_OFFSET..min_len], + §ion_bytes[time_range_offset..min_len], time_range_count, - ), + )?, }) } } @@ -654,19 +739,28 @@ fn parse_reference_time(bytes: &[u8], section: u8) -> Result { fn parse_statistical_time_ranges( bytes: &[u8], time_range_count: usize, -) -> Vec { - bytes - .chunks_exact(12) - .take(time_range_count) - .map(|range| StatisticalTimeRange { +) -> Result> { + let mut ranges = Vec::new(); + ranges + .try_reserve_exact(time_range_count) + .map_err(|error| { + Error::allocation( + "statistical time-range descriptors", + time_range_count, + error, + ) + })?; + for range in bytes.chunks_exact(12).take(time_range_count) { + ranges.push(StatisticalTimeRange { type_of_statistical_processing: range[0], type_of_time_increment: range[1], time_range_unit: range[2], time_range_length: u32::from_be_bytes(range[3..7].try_into().unwrap()), time_increment_unit: range[7], time_increment: u32::from_be_bytes(range[8..12].try_into().unwrap()), - }) - .collect() + }); + } + Ok(ranges) } fn parse_probability_limit(bytes: &[u8]) -> Option { @@ -939,9 +1033,12 @@ mod tests { ); match product.template { ProductDefinitionTemplate::StatisticalProcess(template) => { - assert_eq!(template.time_ranges.len(), 1); - assert_eq!(template.time_ranges[0].type_of_statistical_processing, 1); - assert_eq!(template.time_ranges[0].time_range_length, 6); + assert_eq!(template.interval.time_ranges.len(), 1); + assert_eq!( + template.interval.time_ranges[0].type_of_statistical_processing, + 1 + ); + assert_eq!(template.interval.time_ranges[0].time_range_length, 6); } other => panic!("expected template 4.8, got {other:?}"), } @@ -974,12 +1071,76 @@ mod tests { match product.template { ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { assert_eq!(template.ensemble.perturbation_number, 3); - assert_eq!(template.time_ranges.len(), 1); + assert_eq!(template.interval.time_ranges.len(), 1); } other => panic!("expected template 4.11, got {other:?}"), } } + #[test] + fn parses_probability_statistical_process_template() { + let mut section = product_section_template_zero(); + section.resize(71, 0xff); + set_product_template(&mut section, 9); + section[34] = 1; + section[35] = 10; + section[36] = 3; + section[37] = 1; + section[38..42].copy_from_slice(&crate::binary::encode_wmo_i32(125).unwrap()); + set_statistical_interval(&mut section, 47); + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::ProbabilityStatisticalProcess(template) = product.template + else { + panic!("expected template 4.9"); + }; + assert_eq!( + template.probability.probability, + ProbabilityType::AboveLowerLimit(super::ProbabilityLimit { + scale_factor: 1, + scaled_value: 125, + }) + ); + assert_eq!(template.interval.time_ranges.len(), 1); + assert_eq!(template.interval.end_of_overall_time_interval.hour, 18); + } + + #[test] + fn parses_percentile_statistical_process_template() { + let mut section = product_section_template_zero(); + section.resize(59, 0); + set_product_template(&mut section, 10); + section[34] = 75; + set_statistical_interval(&mut section, 35); + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::PercentileStatisticalProcess(template) = product.template + else { + panic!("expected template 4.10"); + }; + assert_eq!(template.percentile.percentile_value, 75); + assert_eq!(template.interval.time_ranges.len(), 1); + } + + #[test] + fn parses_derived_statistical_process_template() { + let mut section = product_section_template_zero(); + section.resize(60, 0); + set_product_template(&mut section, 12); + section[34] = 1; + section[35] = 20; + set_statistical_interval(&mut section, 36); + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::DerivedStatisticalProcess(template) = product.template + else { + panic!("expected template 4.12"); + }; + assert_eq!(template.derived.derived_forecast_type, 1); + assert_eq!(template.derived.number_of_forecasts_in_ensemble, 20); + assert_eq!(template.interval.time_ranges.len(), 1); + } + #[test] fn rejects_invalid_statistical_process_end_time() { let mut section = product_section_template_eight(); @@ -1052,6 +1213,19 @@ mod tests { } } + #[test] + fn rejects_truncated_interval_product_templates() { + for (template, length) in [(8, 45), (9, 58), (10, 46), (11, 48), (12, 47)] { + let mut section = product_section_template_zero(); + section.resize(length, 0); + set_product_template(&mut section, template); + assert!(matches!( + ProductDefinition::parse(§ion), + Err(Error::InvalidSection { section: 4, .. }) + )); + } + } + fn product_section_template_zero() -> Vec { let mut section = vec![0u8; 34]; section[..4].copy_from_slice(&(34u32).to_be_bytes()); @@ -1072,20 +1246,8 @@ mod tests { fn product_section_template_eight() -> Vec { let mut section = product_section_template_zero(); section.resize(58, 0); - section[..4].copy_from_slice(&(58u32).to_be_bytes()); - section[7..9].copy_from_slice(&8u16.to_be_bytes()); - section[34..36].copy_from_slice(&2026u16.to_be_bytes()); - section[36] = 3; - section[37] = 20; - section[38] = 18; - section[39] = 0; - section[40] = 0; - section[41] = 1; - section[46] = 1; - section[47] = 2; - section[48] = 1; - section[49..53].copy_from_slice(&6u32.to_be_bytes()); - section[53] = 255; + set_product_template(&mut section, 8); + set_statistical_interval(&mut section, 34); section } @@ -1095,6 +1257,24 @@ mod tests { section[7..9].copy_from_slice(&template.to_be_bytes()); } + fn set_statistical_interval(section: &mut [u8], offset: usize) { + section[offset..offset + 2].copy_from_slice(&2026u16.to_be_bytes()); + section[offset + 2] = 3; + section[offset + 3] = 20; + section[offset + 4] = 18; + section[offset + 5] = 0; + section[offset + 6] = 0; + section[offset + 7] = 1; + section[offset + 8..offset + 12].copy_from_slice(&0u32.to_be_bytes()); + let range = offset + 12; + section[range] = 1; + section[range + 1] = 2; + section[range + 2] = 1; + section[range + 3..range + 7].copy_from_slice(&6u32.to_be_bytes()); + section[range + 7] = 255; + section[range + 8..range + 12].copy_from_slice(&0u32.to_be_bytes()); + } + fn valid_identification_section() -> Vec { let mut section = vec![0u8; 21]; section[..4].copy_from_slice(&(21u32).to_be_bytes()); diff --git a/grib-reader/src/lib.rs b/grib-reader/src/lib.rs index dcda18d..fd1195f 100644 --- a/grib-reader/src/lib.rs +++ b/grib-reader/src/lib.rs @@ -58,10 +58,13 @@ pub use parameter::{ LOCAL_PARAMETER_TABLE_CSV_HEADER, }; pub use product::{ - AnalysisOrForecastTemplate, DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, - FixedSurface, Identification, IndividualEnsembleForecastTemplate, PercentileForecastTemplate, - ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ScaledValue, StatisticalProcessTemplate, StatisticalTimeRange, + AnalysisOrForecastTemplate, DerivedForecastTemplate, DerivedStatisticalProcessTemplate, + EnsembleStatisticalProcessTemplate, FixedSurface, Identification, + IndividualEnsembleForecastTemplate, PercentileForecastTemplate, + PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, + ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ScaledValue, StatisticalInterval, StatisticalProcessTemplate, + StatisticalTimeRange, }; use std::io::Read; diff --git a/grib-reader/src/product.rs b/grib-reader/src/product.rs index f8ef688..6de1f16 100644 --- a/grib-reader/src/product.rs +++ b/grib-reader/src/product.rs @@ -1,6 +1,9 @@ pub use grib_core::product::{ - AnalysisOrForecastTemplate, DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, - FixedSurface, Identification, IndividualEnsembleForecastTemplate, PercentileForecastTemplate, - ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ScaledValue, StatisticalProcessTemplate, StatisticalTimeRange, + AnalysisOrForecastTemplate, DerivedForecastTemplate, DerivedStatisticalProcessTemplate, + EnsembleStatisticalProcessTemplate, FixedSurface, Identification, + IndividualEnsembleForecastTemplate, PercentileForecastTemplate, + PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, + ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, ScaledValue, StatisticalInterval, StatisticalProcessTemplate, + StatisticalTimeRange, }; diff --git a/grib-reader/tests/integration.rs b/grib-reader/tests/integration.rs index 3c05919..bdef925 100644 --- a/grib-reader/tests/integration.rs +++ b/grib-reader/tests/integration.rs @@ -425,9 +425,12 @@ fn opens_grib2_statistical_process_product_template() { assert_eq!(valid.hour, 18); match &field.product_definition().unwrap().template { ProductDefinitionTemplate::StatisticalProcess(template) => { - assert_eq!(template.time_ranges.len(), 1); - assert_eq!(template.time_ranges[0].type_of_statistical_processing, 1); - assert_eq!(template.time_ranges[0].time_range_length, 6); + assert_eq!(template.interval.time_ranges.len(), 1); + assert_eq!( + template.interval.time_ranges[0].type_of_statistical_processing, + 1 + ); + assert_eq!(template.interval.time_ranges[0].time_range_length, 6); } other => panic!("expected template 4.8, got {other:?}"), } @@ -451,7 +454,7 @@ fn opens_grib2_ensemble_statistical_process_product_template() { match &field.product_definition().unwrap().template { ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { assert_eq!(template.ensemble.perturbation_number, 3); - assert_eq!(template.time_ranges.len(), 1); + assert_eq!(template.interval.time_ranges.len(), 1); } other => panic!("expected template 4.11, got {other:?}"), } diff --git a/grib-writer/src/lib.rs b/grib-writer/src/lib.rs index ef162f3..f690bb9 100644 --- a/grib-writer/src/lib.rs +++ b/grib-writer/src/lib.rs @@ -15,7 +15,7 @@ use grib_core::{ LambertConformalGrid, LatLonGrid, MercatorGrid, PngPackingParams, PolarStereographicGrid, ProbabilityLimit, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, ProjectedGridCore, ReferenceTime, SimplePackingParams, SpatialDifferencingParams, - StatisticalTimeRange, + StatisticalInterval, StatisticalTimeRange, }; pub use grib_core::grib1::ProductDefinition as Grib1ProductDefinition; @@ -2012,16 +2012,13 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu write_u8_be(out, template.percentile_value) } ProductDefinitionTemplate::StatisticalProcess(template) => { - let range_count = checked_time_range_count(template.time_ranges.len())?; + let range_count = checked_time_range_count(template.interval.time_ranges.len())?; let section_length = statistical_product_section_len(46, range_count)?; write_product_template_prefix(out, product, 8, section_length, &template.base)?; - write_reference_time(out, template.end_of_overall_time_interval)?; - write_u8_be(out, range_count)?; - write_u32_be(out, template.number_of_missing_in_statistical_process)?; - write_statistical_time_ranges(out, &template.time_ranges) + write_statistical_interval(out, &template.interval, range_count) } ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { - let range_count = checked_time_range_count(template.time_ranges.len())?; + let range_count = checked_time_range_count(template.interval.time_ranges.len())?; let section_length = statistical_product_section_len(49, range_count)?; write_product_template_prefix( out, @@ -2031,10 +2028,7 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu &template.ensemble.base, )?; write_ensemble_product_extra(out, &template.ensemble)?; - write_reference_time(out, template.end_of_overall_time_interval)?; - write_u8_be(out, range_count)?; - write_u32_be(out, template.number_of_missing_in_statistical_process)?; - write_statistical_time_ranges(out, &template.time_ranges) + write_statistical_interval(out, &template.interval, range_count) } ProductDefinitionTemplate::Unsupported { number, .. } => { Err(Error::UnsupportedProductTemplate(*number)) @@ -2170,6 +2164,17 @@ fn validate_percentile(percentile: u8) -> Result<()> { } } +fn write_statistical_interval( + out: &mut Vec, + interval: &StatisticalInterval, + range_count: u8, +) -> Result<()> { + write_reference_time(out, interval.end_of_overall_time_interval)?; + write_u8_be(out, range_count)?; + write_u32_be(out, interval.number_of_missing_in_statistical_process)?; + write_statistical_time_ranges(out, &interval.time_ranges) +} + fn write_reference_time(out: &mut Vec, reference_time: ReferenceTime) -> Result<()> { validate_reference_time(reference_time)?; @@ -2526,13 +2531,11 @@ fn validate_supported_product(product: &ProductDefinition) -> Result<()> { } ProductDefinitionTemplate::StatisticalProcess(template) => { validate_product_template_prefix(&template.base)?; - checked_time_range_count(template.time_ranges.len())?; - validate_reference_time(template.end_of_overall_time_interval) + validate_statistical_interval(&template.interval) } ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { validate_product_template_prefix(&template.ensemble.base)?; - checked_time_range_count(template.time_ranges.len())?; - validate_reference_time(template.end_of_overall_time_interval) + validate_statistical_interval(&template.interval) } ProductDefinitionTemplate::Unsupported { number, .. } => { Err(Error::UnsupportedProductTemplate(*number)) @@ -2541,6 +2544,11 @@ fn validate_supported_product(product: &ProductDefinition) -> Result<()> { } } +fn validate_statistical_interval(interval: &StatisticalInterval) -> Result<()> { + checked_time_range_count(interval.time_ranges.len())?; + validate_reference_time(interval.end_of_overall_time_interval) +} + #[cfg(test)] mod tests { use super::{ @@ -2557,8 +2565,8 @@ mod tests { Identification, IndividualEnsembleForecastTemplate, LambertConformalGrid, LatLonGrid, MercatorGrid, PercentileForecastTemplate, PolarStereographicGrid, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ProjectedGridCore, StatisticalProcessTemplate, - StatisticalTimeRange, + ProductDefinitionTemplate, ProjectedGridCore, StatisticalInterval, + StatisticalProcessTemplate, StatisticalTimeRange, }; use grib_reader::sections::scan_sections; use grib_reader::{GribFile, PredefinedBitmap}; @@ -2804,6 +2812,14 @@ mod tests { } } + fn statistical_interval() -> StatisticalInterval { + StatisticalInterval { + end_of_overall_time_interval: interval_end_time(), + number_of_missing_in_statistical_process: 0, + time_ranges: vec![statistical_time_range()], + } + } + fn write_message(fields: impl IntoIterator) -> Vec { let mut bytes = Vec::new(); GribWriter::new(&mut bytes) @@ -3311,9 +3327,7 @@ mod tests { template: ProductDefinitionTemplate::StatisticalProcess( StatisticalProcessTemplate { base, - end_of_overall_time_interval: interval_end_time(), - number_of_missing_in_statistical_process: 0, - time_ranges: vec![statistical_time_range()], + interval: statistical_interval(), }, ), }) @@ -3338,8 +3352,7 @@ mod tests { match &product.template { ProductDefinitionTemplate::StatisticalProcess(template) => { assert_eq!(template.base.forecast_time, 1); - assert_eq!(template.end_of_overall_time_interval, interval_end_time()); - assert_eq!(template.time_ranges, vec![statistical_time_range()]); + assert_eq!(template.interval, statistical_interval()); } other => panic!("expected template 4.8, got {other:?}"), } @@ -3363,9 +3376,7 @@ mod tests { perturbation_number: 3, number_of_forecasts_in_ensemble: 30, }, - end_of_overall_time_interval: interval_end_time(), - number_of_missing_in_statistical_process: 0, - time_ranges: vec![statistical_time_range()], + interval: statistical_interval(), }, ), }) @@ -3392,8 +3403,7 @@ mod tests { assert_eq!(template.ensemble.type_of_ensemble_forecast, 1); assert_eq!(template.ensemble.perturbation_number, 3); assert_eq!(template.ensemble.number_of_forecasts_in_ensemble, 30); - assert_eq!(template.end_of_overall_time_interval, interval_end_time()); - assert_eq!(template.time_ranges, vec![statistical_time_range()]); + assert_eq!(template.interval, statistical_interval()); } other => panic!("expected template 4.11, got {other:?}"), } @@ -3411,9 +3421,10 @@ mod tests { template: ProductDefinitionTemplate::StatisticalProcess( StatisticalProcessTemplate { base: analysis_or_forecast_template(), - end_of_overall_time_interval: interval_end_time(), - number_of_missing_in_statistical_process: 0, - time_ranges: vec![statistical_time_range(); 256], + interval: StatisticalInterval { + time_ranges: vec![statistical_time_range(); 256], + ..statistical_interval() + }, }, ), }) @@ -3438,16 +3449,17 @@ mod tests { template: ProductDefinitionTemplate::StatisticalProcess( StatisticalProcessTemplate { base: analysis_or_forecast_template(), - end_of_overall_time_interval: ReferenceTime { - year: 2026, - month: 13, - day: 20, - hour: 18, - minute: 0, - second: 0, + interval: StatisticalInterval { + end_of_overall_time_interval: ReferenceTime { + year: 2026, + month: 13, + day: 20, + hour: 18, + minute: 0, + second: 0, + }, + ..statistical_interval() }, - number_of_missing_in_statistical_process: 0, - time_ranges: vec![statistical_time_range()], }, ), }) From 5f93aa64048f20fd02ededd9c8091b40075b9ff7 Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 13:53:33 -0500 Subject: [PATCH 6/7] Write interval templates 4.9 4.10 and 4.12 and verify ecCodes metadata --- CHANGELOG.md | 2 + README.md | 4 +- .../fuzz_targets/fuzz_grib_writer_inputs.rs | 68 ++++- grib-writer/src/lib.rs | 186 ++++++++++++- grib-writer/tests/common/mod.rs | 246 +++++++++++++++--- tools/eccodes-reference.c | 155 +++++------ 6 files changed, 509 insertions(+), 152 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 58d87a0..9cde031 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -30,6 +30,8 @@ metadata with ecCodes in the parity suite - decode interval probability, percentile, and derived-ensemble templates 4.9, 4.10, and 4.12 through one shared statistical-interval model +- write interval product templates 4.9, 4.10, and 4.12 and verify their + timestamps, time ranges, thresholds, and ensemble metadata against ecCodes - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index 9b1a554..a6b4fa6 100644 --- a/README.md +++ b/README.md @@ -179,7 +179,7 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Reader GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, and 4.12 - Writer GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, - and 4.11 + 4.9, 4.10, 4.11, and 4.12 - Forecast valid-time helpers for supported fixed-width GRIB1/GRIB2 time units - `GribFile::builder()` for strict or tolerant scanning and allocation limits - Bitmap application with missing values surfaced as `NaN` @@ -204,7 +204,7 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Remaining reader GRIB2 product definition templates beyond 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, and 4.12 - Remaining writer GRIB2 product definition templates beyond 4.0, 4.1, 4.2, - 4.5, 4.6, 4.8, and 4.11 + 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, and 4.12 - Writer GRIB2 row-by-row complex packing Unsupported decode and encode operations fail explicitly with typed errors; diff --git a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs index 1c692af..f8875c1 100644 --- a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs +++ b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs @@ -2,9 +2,11 @@ use grib_core::metadata::ReferenceTime; use grib_core::{ - AnalysisOrForecastTemplate, DerivedForecastTemplate, FixedSurface, GridDefinition, - Identification, LatLonGrid, PercentileForecastTemplate, ProbabilityForecastTemplate, - ProbabilityLimit, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, + AnalysisOrForecastTemplate, DerivedForecastTemplate, DerivedStatisticalProcessTemplate, + FixedSurface, GridDefinition, Identification, LatLonGrid, PercentileForecastTemplate, + PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, + ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, StatisticalInterval, StatisticalTimeRange, }; use grib_reader::GribFile; use grib_writer::{ @@ -225,7 +227,7 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { first_surface: Some(FixedSurface::with_value(103, 0, 850)), second_surface: None, }; - let template = match input.u8() % 4 { + let template = match input.u8() % 7 { 0 => ProductDefinitionTemplate::AnalysisOrForecast(base), 1 => ProductDefinitionTemplate::DerivedForecast(DerivedForecastTemplate { base, @@ -238,10 +240,40 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { total_number_of_forecast_probabilities: input.u8(), probability: probability_type(input), }), - _ => ProductDefinitionTemplate::PercentileForecast(PercentileForecastTemplate { + 3 => ProductDefinitionTemplate::PercentileForecast(PercentileForecastTemplate { base, percentile_value: input.u8() % 101, }), + 4 => ProductDefinitionTemplate::ProbabilityStatisticalProcess( + ProbabilityStatisticalProcessTemplate { + probability: ProbabilityForecastTemplate { + base, + forecast_probability_number: input.u8(), + total_number_of_forecast_probabilities: input.u8(), + probability: probability_type(input), + }, + interval: statistical_interval(input), + }, + ), + 5 => ProductDefinitionTemplate::PercentileStatisticalProcess( + PercentileStatisticalProcessTemplate { + percentile: PercentileForecastTemplate { + base, + percentile_value: input.u8() % 101, + }, + interval: statistical_interval(input), + }, + ), + _ => ProductDefinitionTemplate::DerivedStatisticalProcess( + DerivedStatisticalProcessTemplate { + derived: DerivedForecastTemplate { + base, + derived_forecast_type: input.u8(), + number_of_forecasts_in_ensemble: input.u8(), + }, + interval: statistical_interval(input), + }, + ), }; ProductDefinition { parameter_category, @@ -250,6 +282,32 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { } } +fn statistical_interval(input: &mut Input<'_>) -> StatisticalInterval { + let range_count = usize::from(input.u8() % 4); + let time_ranges = (0..range_count) + .map(|_| StatisticalTimeRange { + type_of_statistical_processing: input.u8(), + type_of_time_increment: input.u8(), + time_range_unit: input.u8(), + time_range_length: u32::from(input.u8()), + time_increment_unit: input.u8(), + time_increment: u32::from(input.u8()), + }) + .collect(); + StatisticalInterval { + end_of_overall_time_interval: ReferenceTime { + year: 2026, + month: 3, + day: 20, + hour: 18, + minute: 0, + second: 0, + }, + number_of_missing_in_statistical_process: u32::from(input.u8()), + time_ranges, + } +} + fn probability_type(input: &mut Input<'_>) -> ProbabilityType { let limit = |input: &mut Input<'_>| ProbabilityLimit { scale_factor: i16::from(input.u8() % 15) - 7, diff --git a/grib-writer/src/lib.rs b/grib-writer/src/lib.rs index f690bb9..4e57d38 100644 --- a/grib-writer/src/lib.rs +++ b/grib-writer/src/lib.rs @@ -1998,8 +1998,7 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu } ProductDefinitionTemplate::DerivedForecast(template) => { write_product_template_prefix(out, product, 2, 36, &template.base)?; - write_u8_be(out, template.derived_forecast_type)?; - write_u8_be(out, template.number_of_forecasts_in_ensemble) + write_derived_product_extra(out, template) } ProductDefinitionTemplate::ProbabilityForecast(template) => { validate_probability_type(template.probability)?; @@ -2017,6 +2016,34 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu write_product_template_prefix(out, product, 8, section_length, &template.base)?; write_statistical_interval(out, &template.interval, range_count) } + ProductDefinitionTemplate::ProbabilityStatisticalProcess(template) => { + validate_probability_type(template.probability.probability)?; + let range_count = checked_time_range_count(template.interval.time_ranges.len())?; + let section_length = statistical_product_section_len(59, range_count)?; + write_product_template_prefix( + out, + product, + 9, + section_length, + &template.probability.base, + )?; + write_probability_product_extra(out, &template.probability)?; + write_statistical_interval(out, &template.interval, range_count) + } + ProductDefinitionTemplate::PercentileStatisticalProcess(template) => { + validate_percentile(template.percentile.percentile_value)?; + let range_count = checked_time_range_count(template.interval.time_ranges.len())?; + let section_length = statistical_product_section_len(47, range_count)?; + write_product_template_prefix( + out, + product, + 10, + section_length, + &template.percentile.base, + )?; + write_u8_be(out, template.percentile.percentile_value)?; + write_statistical_interval(out, &template.interval, range_count) + } ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { let range_count = checked_time_range_count(template.interval.time_ranges.len())?; let section_length = statistical_product_section_len(49, range_count)?; @@ -2030,6 +2057,19 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu write_ensemble_product_extra(out, &template.ensemble)?; write_statistical_interval(out, &template.interval, range_count) } + ProductDefinitionTemplate::DerivedStatisticalProcess(template) => { + let range_count = checked_time_range_count(template.interval.time_ranges.len())?; + let section_length = statistical_product_section_len(48, range_count)?; + write_product_template_prefix( + out, + product, + 12, + section_length, + &template.derived.base, + )?; + write_derived_product_extra(out, &template.derived)?; + write_statistical_interval(out, &template.interval, range_count) + } ProductDefinitionTemplate::Unsupported { number, .. } => { Err(Error::UnsupportedProductTemplate(*number)) } @@ -2096,6 +2136,14 @@ fn write_ensemble_product_extra( write_u8_be(out, template.number_of_forecasts_in_ensemble) } +fn write_derived_product_extra( + out: &mut Vec, + template: &grib_core::DerivedForecastTemplate, +) -> Result<()> { + write_u8_be(out, template.derived_forecast_type)?; + write_u8_be(out, template.number_of_forecasts_in_ensemble) +} + fn write_probability_product_extra( out: &mut Vec, template: &grib_core::ProbabilityForecastTemplate, @@ -2533,10 +2581,24 @@ fn validate_supported_product(product: &ProductDefinition) -> Result<()> { validate_product_template_prefix(&template.base)?; validate_statistical_interval(&template.interval) } + ProductDefinitionTemplate::ProbabilityStatisticalProcess(template) => { + validate_product_template_prefix(&template.probability.base)?; + validate_probability_type(template.probability.probability)?; + validate_statistical_interval(&template.interval) + } + ProductDefinitionTemplate::PercentileStatisticalProcess(template) => { + validate_product_template_prefix(&template.percentile.base)?; + validate_percentile(template.percentile.percentile_value)?; + validate_statistical_interval(&template.interval) + } ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { validate_product_template_prefix(&template.ensemble.base)?; validate_statistical_interval(&template.interval) } + ProductDefinitionTemplate::DerivedStatisticalProcess(template) => { + validate_product_template_prefix(&template.derived.base)?; + validate_statistical_interval(&template.interval) + } ProductDefinitionTemplate::Unsupported { number, .. } => { Err(Error::UnsupportedProductTemplate(*number)) } @@ -2561,12 +2623,13 @@ mod tests { use grib_core::metadata::ReferenceTime; use grib_core::{ AlbersEqualAreaGrid, AnalysisOrForecastTemplate, DataRepresentation, - DerivedForecastTemplate, EnsembleStatisticalProcessTemplate, FixedSurface, GridDefinition, - Identification, IndividualEnsembleForecastTemplate, LambertConformalGrid, LatLonGrid, - MercatorGrid, PercentileForecastTemplate, PolarStereographicGrid, - ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ProjectedGridCore, StatisticalInterval, - StatisticalProcessTemplate, StatisticalTimeRange, + DerivedForecastTemplate, DerivedStatisticalProcessTemplate, + EnsembleStatisticalProcessTemplate, FixedSurface, GridDefinition, Identification, + IndividualEnsembleForecastTemplate, LambertConformalGrid, LatLonGrid, MercatorGrid, + PercentileForecastTemplate, PercentileStatisticalProcessTemplate, PolarStereographicGrid, + ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, + ProbabilityType, ProductDefinition, ProductDefinitionTemplate, ProjectedGridCore, + StatisticalInterval, StatisticalProcessTemplate, StatisticalTimeRange, }; use grib_reader::sections::scan_sections; use grib_reader::{GribFile, PredefinedBitmap}; @@ -3410,6 +3473,113 @@ mod tests { assert_eq!(message.read_flat_data_as_f64().unwrap(), values); } + #[test] + fn writes_probability_statistical_product_template_readable_by_reader() { + let probability = ProbabilityType::AboveLowerLimit(ProbabilityLimit { + scale_factor: 1, + scaled_value: 125, + }); + let field = + field_with_product_template(ProductDefinitionTemplate::ProbabilityStatisticalProcess( + ProbabilityStatisticalProcessTemplate { + probability: ProbabilityForecastTemplate { + base: analysis_or_forecast_template(), + forecast_probability_number: 1, + total_number_of_forecast_probabilities: 10, + probability, + }, + interval: statistical_interval(), + }, + )); + + let bytes = write_message([field]); + let product_section = scan_sections(&bytes) + .unwrap() + .into_iter() + .find(|section| section.number == 4) + .unwrap(); + assert_eq!(product_section.length, 71); + + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::ProbabilityStatisticalProcess(template) = &product.template + else { + panic!("expected template 4.9"); + }; + assert_eq!(template.probability.probability, probability); + assert_eq!(template.interval, statistical_interval()); + assert_eq!(message.valid_time(), Some(interval_end_time())); + } + + #[test] + fn writes_percentile_statistical_product_template_readable_by_reader() { + let field = + field_with_product_template(ProductDefinitionTemplate::PercentileStatisticalProcess( + PercentileStatisticalProcessTemplate { + percentile: PercentileForecastTemplate { + base: analysis_or_forecast_template(), + percentile_value: 75, + }, + interval: statistical_interval(), + }, + )); + + let bytes = write_message([field]); + let product_section = scan_sections(&bytes) + .unwrap() + .into_iter() + .find(|section| section.number == 4) + .unwrap(); + assert_eq!(product_section.length, 59); + + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::PercentileStatisticalProcess(template) = &product.template + else { + panic!("expected template 4.10"); + }; + assert_eq!(template.percentile.percentile_value, 75); + assert_eq!(template.interval, statistical_interval()); + assert_eq!(message.valid_time(), Some(interval_end_time())); + } + + #[test] + fn writes_derived_statistical_product_template_readable_by_reader() { + let field = + field_with_product_template(ProductDefinitionTemplate::DerivedStatisticalProcess( + DerivedStatisticalProcessTemplate { + derived: DerivedForecastTemplate { + base: analysis_or_forecast_template(), + derived_forecast_type: 1, + number_of_forecasts_in_ensemble: 20, + }, + interval: statistical_interval(), + }, + )); + + let bytes = write_message([field]); + let product_section = scan_sections(&bytes) + .unwrap() + .into_iter() + .find(|section| section.number == 4) + .unwrap(); + assert_eq!(product_section.length, 60); + + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::DerivedStatisticalProcess(template) = &product.template + else { + panic!("expected template 4.12"); + }; + assert_eq!(template.derived.derived_forecast_type, 1); + assert_eq!(template.derived.number_of_forecasts_in_ensemble, 20); + assert_eq!(template.interval, statistical_interval()); + assert_eq!(message.valid_time(), Some(interval_end_time())); + } + #[test] fn rejects_too_many_statistical_time_ranges() { let err = Grib2FieldBuilder::new() diff --git a/grib-writer/tests/common/mod.rs b/grib-writer/tests/common/mod.rs index 14bd490..c098955 100644 --- a/grib-writer/tests/common/mod.rs +++ b/grib-writer/tests/common/mod.rs @@ -5,9 +5,11 @@ use std::process::Command; use grib_core::metadata::ReferenceTime; use grib_core::{ - AnalysisOrForecastTemplate, DerivedForecastTemplate, FixedSurface, GridDefinition, - Identification, LatLonGrid, PercentileForecastTemplate, ProbabilityForecastTemplate, - ProbabilityLimit, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, + AnalysisOrForecastTemplate, DerivedForecastTemplate, DerivedStatisticalProcessTemplate, + FixedSurface, GridDefinition, Identification, LatLonGrid, PercentileForecastTemplate, + PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, + ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, + ProductDefinitionTemplate, StatisticalInterval, StatisticalTimeRange, }; use grib_reader::GribFile; use grib_writer::{ @@ -39,6 +41,20 @@ pub struct ReferenceMessage { pub scale_factor_of_upper_limit: Option, pub scaled_value_of_upper_limit: Option, pub percentile_value: Option, + pub interval_end_year: Option, + pub interval_end_month: Option, + pub interval_end_day: Option, + pub interval_end_hour: Option, + pub interval_end_minute: Option, + pub interval_end_second: Option, + pub number_of_time_ranges: Option, + pub number_missing_in_statistical_process: Option, + pub type_of_statistical_processing: Option, + pub type_of_time_increment: Option, + pub time_range_unit: Option, + pub time_range_length: Option, + pub time_increment_unit: Option, + pub time_increment: Option, pub values: Vec>, } @@ -230,49 +246,130 @@ fn assert_product_metadata( match &product.template { ProductDefinitionTemplate::DerivedForecast(template) => { - assert_eq!( - expected.derived_forecast, - Some(i64::from(template.derived_forecast_type)) - ); - assert_eq!( - expected.number_of_forecasts_in_ensemble, - Some(i64::from(template.number_of_forecasts_in_ensemble)) - ); + assert_derived_product_metadata(template, expected) } ProductDefinitionTemplate::ProbabilityForecast(template) => { - assert_eq!( - expected.forecast_probability_number, - Some(i64::from(template.forecast_probability_number)) - ); - assert_eq!( - expected.total_number_of_forecast_probabilities, - Some(i64::from(template.total_number_of_forecast_probabilities)) - ); - assert_eq!( - expected.probability_type, - Some(i64::from(template.probability.code())) - ); - assert_probability_limit_metadata( - template.probability.lower_limit(), - expected.scale_factor_of_lower_limit, - expected.scaled_value_of_lower_limit, - ); - assert_probability_limit_metadata( - template.probability.upper_limit(), - expected.scale_factor_of_upper_limit, - expected.scaled_value_of_upper_limit, - ); + assert_probability_product_metadata(template, expected) } ProductDefinitionTemplate::PercentileForecast(template) => { - assert_eq!( - expected.percentile_value, - Some(i64::from(template.percentile_value)) - ); + assert_percentile_product_metadata(template, expected) + } + ProductDefinitionTemplate::StatisticalProcess(template) => { + assert_statistical_interval_metadata(&template.interval, expected) + } + ProductDefinitionTemplate::ProbabilityStatisticalProcess(template) => { + assert_probability_product_metadata(&template.probability, expected); + assert_statistical_interval_metadata(&template.interval, expected); + } + ProductDefinitionTemplate::PercentileStatisticalProcess(template) => { + assert_percentile_product_metadata(&template.percentile, expected); + assert_statistical_interval_metadata(&template.interval, expected); + } + ProductDefinitionTemplate::EnsembleStatisticalProcess(template) => { + assert_statistical_interval_metadata(&template.interval, expected) + } + ProductDefinitionTemplate::DerivedStatisticalProcess(template) => { + assert_derived_product_metadata(&template.derived, expected); + assert_statistical_interval_metadata(&template.interval, expected); } _ => {} } } +fn assert_derived_product_metadata(actual: &DerivedForecastTemplate, expected: &ReferenceMessage) { + assert_eq!( + expected.derived_forecast, + Some(i64::from(actual.derived_forecast_type)) + ); + assert_eq!( + expected.number_of_forecasts_in_ensemble, + Some(i64::from(actual.number_of_forecasts_in_ensemble)) + ); +} + +fn assert_probability_product_metadata( + actual: &ProbabilityForecastTemplate, + expected: &ReferenceMessage, +) { + assert_eq!( + expected.forecast_probability_number, + Some(i64::from(actual.forecast_probability_number)) + ); + assert_eq!( + expected.total_number_of_forecast_probabilities, + Some(i64::from(actual.total_number_of_forecast_probabilities)) + ); + assert_eq!( + expected.probability_type, + Some(i64::from(actual.probability.code())) + ); + assert_probability_limit_metadata( + actual.probability.lower_limit(), + expected.scale_factor_of_lower_limit, + expected.scaled_value_of_lower_limit, + ); + assert_probability_limit_metadata( + actual.probability.upper_limit(), + expected.scale_factor_of_upper_limit, + expected.scaled_value_of_upper_limit, + ); +} + +fn assert_percentile_product_metadata( + actual: &PercentileForecastTemplate, + expected: &ReferenceMessage, +) { + assert_eq!( + expected.percentile_value, + Some(i64::from(actual.percentile_value)) + ); +} + +fn assert_statistical_interval_metadata(actual: &StatisticalInterval, expected: &ReferenceMessage) { + let end = actual.end_of_overall_time_interval; + assert_eq!(expected.interval_end_year, Some(i64::from(end.year))); + assert_eq!(expected.interval_end_month, Some(i64::from(end.month))); + assert_eq!(expected.interval_end_day, Some(i64::from(end.day))); + assert_eq!(expected.interval_end_hour, Some(i64::from(end.hour))); + assert_eq!(expected.interval_end_minute, Some(i64::from(end.minute))); + assert_eq!(expected.interval_end_second, Some(i64::from(end.second))); + assert_eq!( + expected.number_of_time_ranges, + Some(i64::try_from(actual.time_ranges.len()).unwrap()) + ); + assert_eq!( + expected.number_missing_in_statistical_process, + Some(i64::from(actual.number_of_missing_in_statistical_process)) + ); + + if let [range] = actual.time_ranges.as_slice() { + assert_eq!( + expected.type_of_statistical_processing, + Some(i64::from(range.type_of_statistical_processing)) + ); + assert_eq!( + expected.type_of_time_increment, + Some(i64::from(range.type_of_time_increment)) + ); + assert_eq!( + expected.time_range_unit, + Some(i64::from(range.time_range_unit)) + ); + assert_eq!( + expected.time_range_length, + Some(i64::from(range.time_range_length)) + ); + assert_eq!( + expected.time_increment_unit, + Some(i64::from(range.time_increment_unit)) + ); + assert_eq!( + expected.time_increment, + Some(i64::from(range.time_increment)) + ); + } +} + fn assert_probability_limit_metadata( actual: Option, expected_scale_factor: Option, @@ -382,15 +479,16 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { number_of_forecasts_in_ensemble: 50, }, )); + let probability_type = ProbabilityType::BelowLowerLimit(ProbabilityLimit { + scale_factor: 1, + scaled_value: 2732, + }); let probability = product_field(ProductDefinitionTemplate::ProbabilityForecast( ProbabilityForecastTemplate { base: analysis_or_forecast_template(), forecast_probability_number: 1, total_number_of_forecast_probabilities: 10, - probability: ProbabilityType::BelowLowerLimit(ProbabilityLimit { - scale_factor: 1, - scaled_value: 2732, - }), + probability: probability_type, }, )); let percentile = product_field(ProductDefinitionTemplate::PercentileForecast( @@ -399,6 +497,38 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { percentile_value: 90, }, )); + let probability_interval = + product_field(ProductDefinitionTemplate::ProbabilityStatisticalProcess( + ProbabilityStatisticalProcessTemplate { + probability: ProbabilityForecastTemplate { + base: analysis_or_forecast_template(), + forecast_probability_number: 1, + total_number_of_forecast_probabilities: 10, + probability: probability_type, + }, + interval: statistical_interval(), + }, + )); + let percentile_interval = + product_field(ProductDefinitionTemplate::PercentileStatisticalProcess( + PercentileStatisticalProcessTemplate { + percentile: PercentileForecastTemplate { + base: analysis_or_forecast_template(), + percentile_value: 90, + }, + interval: statistical_interval(), + }, + )); + let derived_interval = product_field(ProductDefinitionTemplate::DerivedStatisticalProcess( + DerivedStatisticalProcessTemplate { + derived: DerivedForecastTemplate { + base: analysis_or_forecast_template(), + derived_forecast_type: 4, + number_of_forecasts_in_ensemble: 50, + }, + interval: statistical_interval(), + }, + )); vec![ ( @@ -420,6 +550,18 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { write_grib2_message([probability]), ), ("writer-percentile.grib2", write_grib2_message([percentile])), + ( + "writer-probability-interval.grib2", + write_grib2_message([probability_interval]), + ), + ( + "writer-percentile-interval.grib2", + write_grib2_message([percentile_interval]), + ), + ( + "writer-derived-interval.grib2", + write_grib2_message([derived_interval]), + ), ("writer-complex.grib2", write_grib2_message([complex])), ( "writer-complex-spatial-first.grib2", @@ -575,6 +717,28 @@ pub fn analysis_or_forecast_template() -> AnalysisOrForecastTemplate { } } +pub fn statistical_interval() -> StatisticalInterval { + StatisticalInterval { + end_of_overall_time_interval: ReferenceTime { + year: 2026, + month: 3, + day: 20, + hour: 18, + minute: 0, + second: 0, + }, + number_of_missing_in_statistical_process: 0, + time_ranges: vec![StatisticalTimeRange { + type_of_statistical_processing: 1, + type_of_time_increment: 2, + time_range_unit: 1, + time_range_length: 6, + time_increment_unit: 255, + time_increment: 0, + }], + } +} + pub fn latlon_grid(ni: u32, nj: u32, scanning_mode: u8) -> GridDefinition { let lon_first = -120_000_000; let lat_first = 50_000_000; diff --git a/tools/eccodes-reference.c b/tools/eccodes-reference.c index 18b3fbf..2c57660 100644 --- a/tools/eccodes-reference.c +++ b/tools/eccodes-reference.c @@ -15,6 +15,13 @@ typedef struct { double checksum; } decode_totals; +typedef struct { + const char *eccodes_key; + const char *json_key; + long value; + int present; +} optional_long_field; + static void die_errno(const char *context, const char *path) { fprintf(stderr, "%s %s: %s\n", context, path, strerror(errno)); exit(1); @@ -76,6 +83,18 @@ static int get_optional_long( return 1; } +static void load_optional_long_fields( + codes_handle *handle, + optional_long_field *fields, + size_t field_count, + const char *path +) { + for (size_t i = 0; i < field_count; ++i) { + fields[i].present = + get_optional_long(handle, fields[i].eccodes_key, &fields[i].value, path); + } +} + static long get_grid_dimension( codes_handle *handle, const char *primary_key, @@ -278,6 +297,15 @@ static void print_optional_long(int present, long value) { } } +static void print_optional_long_fields(const optional_long_field *fields, size_t field_count) { + for (size_t i = 0; i < field_count; ++i) { + fputc(',', stdout); + print_json_string(fields[i].json_key); + fputc(':', stdout); + print_optional_long(fields[i].present, fields[i].value); + } +} + static decode_totals decode_file(const char *path, int emit_json) { FILE *fp = fopen(path, "rb"); if (fp == NULL) { @@ -309,71 +337,36 @@ static decode_totals decode_file(const char *path, int emit_json) { long second = get_long(handle, "second", path); long ni = get_grid_dimension(handle, "Ni", "Nx", path); long nj = get_grid_dimension(handle, "Nj", "Ny", path); - long product_definition_template_number = 0; - long derived_forecast = 0; - long number_of_forecasts_in_ensemble = 0; - long forecast_probability_number = 0; - long total_number_of_forecast_probabilities = 0; - long probability_type = 0; - long scale_factor_of_lower_limit = 0; - long scaled_value_of_lower_limit = 0; - long scale_factor_of_upper_limit = 0; - long scaled_value_of_upper_limit = 0; - long percentile_value = 0; - int has_product_definition_template_number = get_optional_long( - handle, - "productDefinitionTemplateNumber", - &product_definition_template_number, - path - ); - int has_derived_forecast = - get_optional_long(handle, "derivedForecast", &derived_forecast, path); - int has_number_of_forecasts_in_ensemble = get_optional_long( - handle, - "numberOfForecastsInEnsemble", - &number_of_forecasts_in_ensemble, - path - ); - int has_forecast_probability_number = get_optional_long( - handle, - "forecastProbabilityNumber", - &forecast_probability_number, - path - ); - int has_total_number_of_forecast_probabilities = get_optional_long( - handle, - "totalNumberOfForecastProbabilities", - &total_number_of_forecast_probabilities, - path - ); - int has_probability_type = - get_optional_long(handle, "probabilityType", &probability_type, path); - int has_scale_factor_of_lower_limit = get_optional_long( - handle, - "scaleFactorOfLowerLimit", - &scale_factor_of_lower_limit, - path - ); - int has_scaled_value_of_lower_limit = get_optional_long( - handle, - "scaledValueOfLowerLimit", - &scaled_value_of_lower_limit, - path - ); - int has_scale_factor_of_upper_limit = get_optional_long( - handle, - "scaleFactorOfUpperLimit", - &scale_factor_of_upper_limit, - path - ); - int has_scaled_value_of_upper_limit = get_optional_long( - handle, - "scaledValueOfUpperLimit", - &scaled_value_of_upper_limit, - path - ); - int has_percentile_value = - get_optional_long(handle, "percentileValue", &percentile_value, path); + optional_long_field product_metadata[] = { + {"productDefinitionTemplateNumber", "product_definition_template_number", 0, 0}, + {"derivedForecast", "derived_forecast", 0, 0}, + {"numberOfForecastsInEnsemble", "number_of_forecasts_in_ensemble", 0, 0}, + {"forecastProbabilityNumber", "forecast_probability_number", 0, 0}, + {"totalNumberOfForecastProbabilities", "total_number_of_forecast_probabilities", 0, 0}, + {"probabilityType", "probability_type", 0, 0}, + {"scaleFactorOfLowerLimit", "scale_factor_of_lower_limit", 0, 0}, + {"scaledValueOfLowerLimit", "scaled_value_of_lower_limit", 0, 0}, + {"scaleFactorOfUpperLimit", "scale_factor_of_upper_limit", 0, 0}, + {"scaledValueOfUpperLimit", "scaled_value_of_upper_limit", 0, 0}, + {"percentileValue", "percentile_value", 0, 0}, + {"yearOfEndOfOverallTimeInterval", "interval_end_year", 0, 0}, + {"monthOfEndOfOverallTimeInterval", "interval_end_month", 0, 0}, + {"dayOfEndOfOverallTimeInterval", "interval_end_day", 0, 0}, + {"hourOfEndOfOverallTimeInterval", "interval_end_hour", 0, 0}, + {"minuteOfEndOfOverallTimeInterval", "interval_end_minute", 0, 0}, + {"secondOfEndOfOverallTimeInterval", "interval_end_second", 0, 0}, + {"numberOfTimeRange", "number_of_time_ranges", 0, 0}, + {"numberOfMissingInStatisticalProcess", "number_missing_in_statistical_process", 0, 0}, + {"typeOfStatisticalProcessing", "type_of_statistical_processing", 0, 0}, + {"typeOfTimeIncrement", "type_of_time_increment", 0, 0}, + {"indicatorOfUnitForTimeRange", "time_range_unit", 0, 0}, + {"lengthOfTimeRange", "time_range_length", 0, 0}, + {"indicatorOfUnitForTimeIncrement", "time_increment_unit", 0, 0}, + {"timeIncrement", "time_increment", 0, 0}, + }; + const size_t product_metadata_count = + sizeof(product_metadata) / sizeof(product_metadata[0]); + load_optional_long_fields(handle, product_metadata, product_metadata_count, path); char name[256]; get_string(handle, "name", name, sizeof(name), path); @@ -412,37 +405,7 @@ static decode_totals decode_file(const char *path, int emit_json) { fprintf(stdout, "%ld", ni); fputs(",\"nj\":", stdout); fprintf(stdout, "%ld", nj); - fputs(",\"product_definition_template_number\":", stdout); - print_optional_long( - has_product_definition_template_number, - product_definition_template_number - ); - fputs(",\"derived_forecast\":", stdout); - print_optional_long(has_derived_forecast, derived_forecast); - fputs(",\"number_of_forecasts_in_ensemble\":", stdout); - print_optional_long( - has_number_of_forecasts_in_ensemble, - number_of_forecasts_in_ensemble - ); - fputs(",\"forecast_probability_number\":", stdout); - print_optional_long(has_forecast_probability_number, forecast_probability_number); - fputs(",\"total_number_of_forecast_probabilities\":", stdout); - print_optional_long( - has_total_number_of_forecast_probabilities, - total_number_of_forecast_probabilities - ); - fputs(",\"probability_type\":", stdout); - print_optional_long(has_probability_type, probability_type); - fputs(",\"scale_factor_of_lower_limit\":", stdout); - print_optional_long(has_scale_factor_of_lower_limit, scale_factor_of_lower_limit); - fputs(",\"scaled_value_of_lower_limit\":", stdout); - print_optional_long(has_scaled_value_of_lower_limit, scaled_value_of_lower_limit); - fputs(",\"scale_factor_of_upper_limit\":", stdout); - print_optional_long(has_scale_factor_of_upper_limit, scale_factor_of_upper_limit); - fputs(",\"scaled_value_of_upper_limit\":", stdout); - print_optional_long(has_scaled_value_of_upper_limit, scaled_value_of_upper_limit); - fputs(",\"percentile_value\":", stdout); - print_optional_long(has_percentile_value, percentile_value); + print_optional_long_fields(product_metadata, product_metadata_count); fputs(",\"values\":[", stdout); for (size_t i = 0; i < value_len; ++i) { if (i > 0) { From ef9bb4d1acd466d1befb755942994cd86712ed0f Mon Sep 17 00:00:00 2001 From: i-norden Date: Wed, 15 Jul 2026 15:18:55 -0500 Subject: [PATCH 7/7] Decode and write spatial product template 4.15 with ecCodes parity --- CHANGELOG.md | 2 + README.md | 8 +-- grib-core/src/lib.rs | 4 +- grib-core/src/product.rs | 54 ++++++++++++++++++- .../fuzz_targets/fuzz_grib_writer_inputs.rs | 12 +++-- grib-reader/src/lib.rs | 4 +- grib-reader/src/product.rs | 4 +- grib-writer/src/lib.rs | 46 +++++++++++++++- grib-writer/tests/common/mod.rs | 28 +++++++++- tools/eccodes-reference.c | 3 ++ 10 files changed, 148 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9cde031..684eeeb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,8 @@ 4.10, and 4.12 through one shared statistical-interval model - write interval product templates 4.9, 4.10, and 4.12 and verify their timestamps, time ranges, thresholds, and ensemble metadata against ecCodes +- decode and write spatially processed product template 4.15 and verify its + process type, spatial method, and source-point count against ecCodes - raise the workspace MSRV to Rust 1.87 for the coordinated breaking release ## 0.6.0 - 2026-06-25 diff --git a/README.md b/README.md index a6b4fa6..0d6eb30 100644 --- a/README.md +++ b/README.md @@ -177,9 +177,9 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Explicit north-up decode methods in addition to reader order that preserves the encoded i/j scan directions - Reader GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, - 4.9, 4.10, 4.11, and 4.12 + 4.9, 4.10, 4.11, 4.12, and 4.15 - Writer GRIB2 product definition templates 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, - 4.9, 4.10, 4.11, and 4.12 + 4.9, 4.10, 4.11, 4.12, and 4.15 - Forecast valid-time helpers for supported fixed-width GRIB1/GRIB2 time units - `GribFile::builder()` for strict or tolerant scanning and allocation limits - Bitmap application with missing values surfaced as `NaN` @@ -202,9 +202,9 @@ GribWriter::new(&mut bytes).write_grib2_message([field])?; - Remaining GRIB2 grid templates beyond 3.0, 3.1, 3.10, 3.20, 3.30, 3.31, and 3.40 - Remaining reader GRIB2 product definition templates beyond 4.0, 4.1, 4.2, - 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, and 4.12 + 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, 4.12, and 4.15 - Remaining writer GRIB2 product definition templates beyond 4.0, 4.1, 4.2, - 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, and 4.12 + 4.5, 4.6, 4.8, 4.9, 4.10, 4.11, 4.12, and 4.15 - Writer GRIB2 row-by-row complex packing Unsupported decode and encode operations fail explicitly with typed errors; diff --git a/grib-core/src/lib.rs b/grib-core/src/lib.rs index 0ce23c5..7052baa 100644 --- a/grib-core/src/lib.rs +++ b/grib-core/src/lib.rs @@ -34,6 +34,6 @@ pub use product::{ IndividualEnsembleForecastTemplate, PercentileForecastTemplate, PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ScaledValue, StatisticalInterval, StatisticalProcessTemplate, - StatisticalTimeRange, + ProductDefinitionTemplate, ScaledValue, SpatialProcessTemplate, StatisticalInterval, + StatisticalProcessTemplate, StatisticalTimeRange, }; diff --git a/grib-core/src/product.rs b/grib-core/src/product.rs index 750763e..bd3c5b6 100644 --- a/grib-core/src/product.rs +++ b/grib-core/src/product.rs @@ -129,6 +129,7 @@ pub enum ProductDefinitionTemplate { PercentileStatisticalProcess(PercentileStatisticalProcessTemplate), EnsembleStatisticalProcess(EnsembleStatisticalProcessTemplate), DerivedStatisticalProcess(DerivedStatisticalProcessTemplate), + SpatialProcess(SpatialProcessTemplate), /// A well-framed Section 4 whose template is not interpreted by this /// version of the library. `raw` contains the template-specific bytes /// following the common parameter category and number. @@ -334,6 +335,15 @@ pub struct DerivedStatisticalProcessTemplate { pub interval: StatisticalInterval, } +/// Product Definition Template 4.15: statistical processing over a spatial area. +#[derive(Debug, Clone, PartialEq)] +pub struct SpatialProcessTemplate { + pub base: AnalysisOrForecastTemplate, + pub statistical_process: u8, + pub spatial_processing: u8, + pub number_of_points_used: u8, +} + /// Statistical processing descriptor from GRIB2 Product Definition templates /// with one or more time range specifications. #[derive(Debug, Clone, PartialEq, Eq)] @@ -454,6 +464,9 @@ impl ProductDefinitionTemplate { 12 => Ok(Self::DerivedStatisticalProcess( DerivedStatisticalProcessTemplate::parse(section_bytes)?, )), + 15 => Ok(Self::SpatialProcess(SpatialProcessTemplate::parse( + section_bytes, + )?)), number => { let raw_len = section_bytes.len() - 11; let mut raw = Vec::new(); @@ -478,6 +491,7 @@ impl ProductDefinitionTemplate { Self::PercentileStatisticalProcess(_) => 10, Self::EnsembleStatisticalProcess(_) => 11, Self::DerivedStatisticalProcess(_) => 12, + Self::SpatialProcess(_) => 15, Self::Unsupported { number, .. } => *number, } } @@ -494,6 +508,7 @@ impl ProductDefinitionTemplate { Self::PercentileStatisticalProcess(template) => &template.percentile.base, Self::EnsembleStatisticalProcess(template) => &template.ensemble.base, Self::DerivedStatisticalProcess(template) => &template.derived.base, + Self::SpatialProcess(template) => &template.base, Self::Unsupported { .. } => return None, }) } @@ -519,7 +534,8 @@ impl ProductDefinitionTemplate { | Self::IndividualEnsembleForecast(_) | Self::DerivedForecast(_) | Self::ProbabilityForecast(_) - | Self::PercentileForecast(_) => None, + | Self::PercentileForecast(_) + | Self::SpatialProcess(_) => None, Self::Unsupported { .. } => None, } } @@ -665,6 +681,21 @@ impl DerivedStatisticalProcessTemplate { } } +impl SpatialProcessTemplate { + const MINIMUM_LENGTH: usize = 37; + + fn parse(section_bytes: &[u8]) -> Result { + require_len(section_bytes, Self::MINIMUM_LENGTH, "template 4.15")?; + + Ok(Self { + base: AnalysisOrForecastTemplate::parse(section_bytes)?, + statistical_process: section_bytes[34], + spatial_processing: section_bytes[35], + number_of_points_used: section_bytes[36], + }) + } +} + impl StatisticalInterval { fn parse(section_bytes: &[u8], end_time_offset: usize, context: &str) -> Result { let time_range_offset = @@ -1141,6 +1172,25 @@ mod tests { assert_eq!(template.interval.time_ranges.len(), 1); } + #[test] + fn parses_spatial_process_template() { + let mut section = product_section_template_zero(); + section.resize(37, 0); + set_product_template(&mut section, 15); + section[34] = 1; + section[35] = 2; + section[36] = 4; + + let product = ProductDefinition::parse(§ion).unwrap(); + let ProductDefinitionTemplate::SpatialProcess(template) = product.template else { + panic!("expected template 4.15"); + }; + assert_eq!(template.statistical_process, 1); + assert_eq!(template.spatial_processing, 2); + assert_eq!(template.number_of_points_used, 4); + assert_eq!(template.base.forecast_time, 6); + } + #[test] fn rejects_invalid_statistical_process_end_time() { let mut section = product_section_template_eight(); @@ -1202,7 +1252,7 @@ mod tests { #[test] fn rejects_truncated_instantaneous_product_templates() { - for (template, length) in [(2, 35), (5, 46), (6, 34)] { + for (template, length) in [(2, 35), (5, 46), (6, 34), (15, 36)] { let mut section = product_section_template_zero(); section.resize(length, 0); set_product_template(&mut section, template); diff --git a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs index f8875c1..67a631a 100644 --- a/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs +++ b/grib-reader/fuzz/fuzz_targets/fuzz_grib_writer_inputs.rs @@ -6,7 +6,7 @@ use grib_core::{ FixedSurface, GridDefinition, Identification, LatLonGrid, PercentileForecastTemplate, PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, StatisticalInterval, StatisticalTimeRange, + ProductDefinitionTemplate, SpatialProcessTemplate, StatisticalInterval, StatisticalTimeRange, }; use grib_reader::GribFile; use grib_writer::{ @@ -227,7 +227,7 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { first_surface: Some(FixedSurface::with_value(103, 0, 850)), second_surface: None, }; - let template = match input.u8() % 7 { + let template = match input.u8() % 8 { 0 => ProductDefinitionTemplate::AnalysisOrForecast(base), 1 => ProductDefinitionTemplate::DerivedForecast(DerivedForecastTemplate { base, @@ -264,7 +264,7 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { interval: statistical_interval(input), }, ), - _ => ProductDefinitionTemplate::DerivedStatisticalProcess( + 6 => ProductDefinitionTemplate::DerivedStatisticalProcess( DerivedStatisticalProcessTemplate { derived: DerivedForecastTemplate { base, @@ -274,6 +274,12 @@ fn product(input: &mut Input<'_>) -> ProductDefinition { interval: statistical_interval(input), }, ), + _ => ProductDefinitionTemplate::SpatialProcess(SpatialProcessTemplate { + base, + statistical_process: input.u8(), + spatial_processing: input.u8(), + number_of_points_used: input.u8(), + }), }; ProductDefinition { parameter_category, diff --git a/grib-reader/src/lib.rs b/grib-reader/src/lib.rs index fd1195f..8d55d7d 100644 --- a/grib-reader/src/lib.rs +++ b/grib-reader/src/lib.rs @@ -63,8 +63,8 @@ pub use product::{ IndividualEnsembleForecastTemplate, PercentileForecastTemplate, PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ScaledValue, StatisticalInterval, StatisticalProcessTemplate, - StatisticalTimeRange, + ProductDefinitionTemplate, ScaledValue, SpatialProcessTemplate, StatisticalInterval, + StatisticalProcessTemplate, StatisticalTimeRange, }; use std::io::Read; diff --git a/grib-reader/src/product.rs b/grib-reader/src/product.rs index 6de1f16..9c10647 100644 --- a/grib-reader/src/product.rs +++ b/grib-reader/src/product.rs @@ -4,6 +4,6 @@ pub use grib_core::product::{ IndividualEnsembleForecastTemplate, PercentileForecastTemplate, PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, ScaledValue, StatisticalInterval, StatisticalProcessTemplate, - StatisticalTimeRange, + ProductDefinitionTemplate, ScaledValue, SpatialProcessTemplate, StatisticalInterval, + StatisticalProcessTemplate, StatisticalTimeRange, }; diff --git a/grib-writer/src/lib.rs b/grib-writer/src/lib.rs index 4e57d38..42daf78 100644 --- a/grib-writer/src/lib.rs +++ b/grib-writer/src/lib.rs @@ -2070,6 +2070,12 @@ fn write_product_section(out: &mut Vec, product: &ProductDefinition) -> Resu write_derived_product_extra(out, &template.derived)?; write_statistical_interval(out, &template.interval, range_count) } + ProductDefinitionTemplate::SpatialProcess(template) => { + write_product_template_prefix(out, product, 15, 37, &template.base)?; + write_u8_be(out, template.statistical_process)?; + write_u8_be(out, template.spatial_processing)?; + write_u8_be(out, template.number_of_points_used) + } ProductDefinitionTemplate::Unsupported { number, .. } => { Err(Error::UnsupportedProductTemplate(*number)) } @@ -2599,6 +2605,9 @@ fn validate_supported_product(product: &ProductDefinition) -> Result<()> { validate_product_template_prefix(&template.derived.base)?; validate_statistical_interval(&template.interval) } + ProductDefinitionTemplate::SpatialProcess(template) => { + validate_product_template_prefix(&template.base) + } ProductDefinitionTemplate::Unsupported { number, .. } => { Err(Error::UnsupportedProductTemplate(*number)) } @@ -2629,7 +2638,8 @@ mod tests { PercentileForecastTemplate, PercentileStatisticalProcessTemplate, PolarStereographicGrid, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, ProductDefinitionTemplate, ProjectedGridCore, - StatisticalInterval, StatisticalProcessTemplate, StatisticalTimeRange, + SpatialProcessTemplate, StatisticalInterval, StatisticalProcessTemplate, + StatisticalTimeRange, }; use grib_reader::sections::scan_sections; use grib_reader::{GribFile, PredefinedBitmap}; @@ -3333,6 +3343,40 @@ mod tests { ); } + #[test] + fn writes_spatial_process_product_template_readable_by_reader() { + let field = field_with_product_template(ProductDefinitionTemplate::SpatialProcess( + SpatialProcessTemplate { + base: analysis_or_forecast_template(), + statistical_process: 1, + spatial_processing: 2, + number_of_points_used: 4, + }, + )); + + let bytes = write_message([field]); + let product_section = scan_sections(&bytes) + .unwrap() + .into_iter() + .find(|section| section.number == 4) + .unwrap(); + assert_eq!(product_section.length, 37); + + let file = GribFile::from_bytes(bytes).unwrap(); + let message = file.message(0).unwrap(); + let product = message.product_definition().unwrap(); + let ProductDefinitionTemplate::SpatialProcess(template) = &product.template else { + panic!("expected template 4.15"); + }; + assert_eq!(template.statistical_process, 1); + assert_eq!(template.spatial_processing, 2); + assert_eq!(template.number_of_points_used, 4); + assert_eq!( + message.read_flat_data_as_f64().unwrap(), + [1.0, 2.0, 3.0, 4.0] + ); + } + #[test] fn rejects_noncanonical_probability_types_and_invalid_percentiles() { for template in [ diff --git a/grib-writer/tests/common/mod.rs b/grib-writer/tests/common/mod.rs index c098955..a376154 100644 --- a/grib-writer/tests/common/mod.rs +++ b/grib-writer/tests/common/mod.rs @@ -9,7 +9,7 @@ use grib_core::{ FixedSurface, GridDefinition, Identification, LatLonGrid, PercentileForecastTemplate, PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit, ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition, - ProductDefinitionTemplate, StatisticalInterval, StatisticalTimeRange, + ProductDefinitionTemplate, SpatialProcessTemplate, StatisticalInterval, StatisticalTimeRange, }; use grib_reader::GribFile; use grib_writer::{ @@ -55,6 +55,9 @@ pub struct ReferenceMessage { pub time_range_length: Option, pub time_increment_unit: Option, pub time_increment: Option, + pub spatial_statistical_process: Option, + pub spatial_processing: Option, + pub number_of_points_used: Option, pub values: Vec>, } @@ -272,6 +275,20 @@ fn assert_product_metadata( assert_derived_product_metadata(&template.derived, expected); assert_statistical_interval_metadata(&template.interval, expected); } + ProductDefinitionTemplate::SpatialProcess(template) => { + assert_eq!( + expected.spatial_statistical_process, + Some(i64::from(template.statistical_process)) + ); + assert_eq!( + expected.spatial_processing, + Some(i64::from(template.spatial_processing)) + ); + assert_eq!( + expected.number_of_points_used, + Some(i64::from(template.number_of_points_used)) + ); + } _ => {} } } @@ -529,6 +546,14 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { interval: statistical_interval(), }, )); + let spatial = product_field(ProductDefinitionTemplate::SpatialProcess( + SpatialProcessTemplate { + base: analysis_or_forecast_template(), + statistical_process: 1, + spatial_processing: 2, + number_of_points_used: 4, + }, + )); vec![ ( @@ -562,6 +587,7 @@ pub fn writer_reference_samples() -> Vec<(&'static str, Vec)> { "writer-derived-interval.grib2", write_grib2_message([derived_interval]), ), + ("writer-spatial.grib2", write_grib2_message([spatial])), ("writer-complex.grib2", write_grib2_message([complex])), ( "writer-complex-spatial-first.grib2", diff --git a/tools/eccodes-reference.c b/tools/eccodes-reference.c index 2c57660..6714843 100644 --- a/tools/eccodes-reference.c +++ b/tools/eccodes-reference.c @@ -363,6 +363,9 @@ static decode_totals decode_file(const char *path, int emit_json) { {"lengthOfTimeRange", "time_range_length", 0, 0}, {"indicatorOfUnitForTimeIncrement", "time_increment_unit", 0, 0}, {"timeIncrement", "time_increment", 0, 0}, + {"statisticalProcess", "spatial_statistical_process", 0, 0}, + {"spatialProcessing", "spatial_processing", 0, 0}, + {"numberOfPointsUsed", "number_of_points_used", 0, 0}, }; const size_t product_metadata_count = sizeof(product_metadata) / sizeof(product_metadata[0]);