Add ensemble wave period template and probability units

gribberish/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "gribberish"
-version = "1.7.0"
+version = "1.8.0"
 authors = ["Matthew Iannucci <mpiannucci@gmail.com>"]
 description = "Parse grib 2 files with Rust"
 edition = "2021"
@@ -11,8 +11,8 @@ categories = ["science", "encoding", "compression"]
 exclude = ["tests"]
 
 [dependencies]
-gribberish-types = { path = "../types", version = "1.7.0" }
-gribberish-macros = { path = "../macros", version = "1.7.0" }
+gribberish-types = { path = "../types", version = "1.8.0" }
+gribberish-macros = { path = "../macros", version = "1.8.0" }
 chrono = "0.4"
 openjpeg-sys = { version = "1.0.3", optional = true }
 png = { version = "0.17.2", optional = true }

gribberish/src/message.rs

@@ -235,7 +235,7 @@ impl<'a> Message<'a> {
             ""
         };
 
-        // Wave period band (template 4.103) - distinguishes otherwise identical
+        // Wave period band (templates 4.103/4.104) - distinguishes otherwise identical
         // period-banded significant wave height messages from one another.
         let wave_period = match self.wave_period_range().unwrap_or(None) {
             Some((lower, upper)) => {
@@ -502,6 +502,12 @@ impl<'a> Message<'a> {
                 Ok(unit)
             }
             Message::Grib2 { .. } => {
+                // Probability products (templates 4.5 and 4.9) carry the
+                // parameter the threshold applies to, but their values are
+                // percentages rather than that parameter's unit.
+                if self.probability_type()?.is_some() {
+                    return Ok("%".to_string());
+                }
                 let parameter = self.parameter()?;
                 Ok(parameter.unit)
             }
@@ -602,7 +608,7 @@ impl<'a> Message<'a> {
     }
 
     /// Returns the inclusive wave period range `(lower, upper)` in seconds for
-    /// messages that select waves by period band (template 4.103). Returns
+    /// messages that select waves by period band (templates 4.103 and 4.104). Returns
     /// `None` for any other product template.
     pub fn wave_period_range(&self) -> Result<Option<WavePeriodRange>, GribberishError> {
         match self {

gribberish/src/sections/product_definition.rs

@@ -9,7 +9,7 @@ use crate::{
         HorizontalAnalysisForecastTemplate, HorizontalEnsembleForecastTemplate,
         PercentileHorizontalTemplate, PercentileHorizontalTimeIntervalTemplate,
         ProbabilityHorizontalForecastTemplate, ProbabilityHorizontalTimeIntervalTemplate,
-        WavePeriodRangeHorizontalForecastTemplate,
+        WavePeriodRangeEnsembleForecastTemplate, WavePeriodRangeHorizontalForecastTemplate,
     },
     utils::{read_u16_from_bytes, read_u32_from_bytes},
 };
@@ -83,6 +83,10 @@ impl<'a> ProductDefinitionSection<'a> {
                 self.data.to_vec(),
                 discipline,
             ))),
+            104 => Some(Box::new(WavePeriodRangeEnsembleForecastTemplate::new(
+                self.data.to_vec(),
+                discipline,
+            ))),
             107 => Some(Box::new(
                 DerivedEnsembleForecastTimeIntervalReferenceTemplate::new(
                     self.data.to_vec(),

gribberish/src/templates/product/mod.rs

@@ -12,6 +12,7 @@ pub mod probability_horizontal_template;
 pub mod probability_horizontal_time_interval_template;
 pub mod product_template;
 pub mod tables;
+pub mod wave_period_range_ensemble_template;
 pub mod wave_period_range_horizontal_template;
 
 pub use average_accumulation_extreme_horizontal_analysis_template::AverageAccumulationExtremeHorizontalAnalysisForecastTemplate;
@@ -24,4 +25,5 @@ pub use percentile_horizontal_template::PercentileHorizontalTemplate;
 pub use percentile_horizontal_time_interval_template::PercentileHorizontalTimeIntervalTemplate;
 pub use probability_horizontal_template::ProbabilityHorizontalForecastTemplate;
 pub use probability_horizontal_time_interval_template::ProbabilityHorizontalTimeIntervalTemplate;
+pub use wave_period_range_ensemble_template::WavePeriodRangeEnsembleForecastTemplate;
 pub use wave_period_range_horizontal_template::WavePeriodRangeHorizontalForecastTemplate;

gribberish/src/templates/product/product_template.rs

@@ -73,7 +73,7 @@ pub trait ProductTemplate {
     }
 
     /// Returns the inclusive wave period range `(lower, upper)` in seconds for
-    /// templates that select waves by period band (e.g. template 4.103). Either
+    /// templates that select waves by period band (templates 4.103 and 4.104). Either
     /// limit may be `None` if the range is open ended. Returns `None` for
     /// templates that do not define a wave period range.
     fn wave_period_range(&self) -> Option<WavePeriodRange> {

gribberish/src/templates/product/wave_period_range_ensemble_template.rs

@@ -0,0 +1,141 @@
+use crate::templates::template::{Template, TemplateType};
+use chrono::{DateTime, Utc};
+
+use super::product_template::{ProductTemplate, WavePeriodRange};
+use super::tables::{EnsembleForecastType, FixedSurfaceType, GeneratingProcess, TimeUnit};
+use super::WavePeriodRangeHorizontalForecastTemplate;
+
+/// GRIB2 Product Definition Template 4.104
+///
+/// "Individual ensemble forecast, control and perturbed, at a horizontal level
+/// or in a horizontal layer at a point in time for waves selected by period
+/// range."
+///
+/// This is template 4.103 ([`WavePeriodRangeHorizontalForecastTemplate`]) with
+/// the three template 4.1 ensemble octets (type of ensemble forecast,
+/// perturbation number, number of forecasts in ensemble) appended after the
+/// second fixed surface. Octets 10-45 share the 4.103 layout, so those fields
+/// are read through a wrapped 4.103 template.
+///
+/// Used by the ECMWF IFS wave model for the period-banded significant wave
+/// heights (`h1012`, `h1214`, ...) in both HRES (perturbation number 0) and ENS.
+pub struct WavePeriodRangeEnsembleForecastTemplate {
+    inner: WavePeriodRangeHorizontalForecastTemplate,
+}
+
+impl Template for WavePeriodRangeEnsembleForecastTemplate {
+    fn data(&self) -> &[u8] {
+        self.inner.data()
+    }
+
+    fn template_number(&self) -> u16 {
+        104
+    }
+
+    fn template_type(&self) -> TemplateType {
+        TemplateType::Product
+    }
+
+    fn template_name(&self) -> &str {
+        "Individual ensemble forecast, control and perturbed, at a horizontal level or in a horizontal layer at a point in time for waves selected by period range"
+    }
+}
+
+impl WavePeriodRangeEnsembleForecastTemplate {
+    pub fn new(data: Vec<u8>, discipline: u8) -> Self {
+        WavePeriodRangeEnsembleForecastTemplate {
+            inner: WavePeriodRangeHorizontalForecastTemplate::new(data, discipline),
+        }
+    }
+
+    /// Period range fields shared with template 4.103.
+    pub fn wave_period_range_template(&self) -> &WavePeriodRangeHorizontalForecastTemplate {
+        &self.inner
+    }
+
+    pub fn type_of_ensemble_forecast(&self) -> EnsembleForecastType {
+        self.data()[45].into()
+    }
+
+    pub fn perturbation_number(&self) -> u8 {
+        self.data()[46]
+    }
+
+    pub fn number_of_forecasts_in_ensemble(&self) -> u8 {
+        self.data()[47]
+    }
+}
+
+impl ProductTemplate for WavePeriodRangeEnsembleForecastTemplate {
+    fn discipline(&self) -> u8 {
+        self.inner.discipline()
+    }
+
+    fn category_value(&self) -> u8 {
+        self.inner.category_value()
+    }
+
+    fn parameter_value(&self) -> u8 {
+        self.inner.parameter_value()
+    }
+
+    fn generating_process(&self) -> GeneratingProcess {
+        self.inner.generating_process()
+    }
+
+    fn time_unit(&self) -> TimeUnit {
+        self.inner.time_unit()
+    }
+
+    fn time_increment_unit(&self) -> Option<TimeUnit> {
+        None
+    }
+
+    fn time_interval(&self) -> i32 {
+        self.inner.time_interval()
+    }
+
+    fn time_increment_interval(&self) -> Option<u32> {
+        None
+    }
+
+    fn forecast_end_datetime(&self, _reference_date: DateTime<Utc>) -> Option<DateTime<Utc>> {
+        None
+    }
+
+    fn first_fixed_surface_type(&self) -> FixedSurfaceType {
+        self.inner.first_fixed_surface_type()
+    }
+
+    fn first_fixed_surface_value(&self) -> Option<f64> {
+        self.inner.first_fixed_surface_value()
+    }
+
+    fn second_fixed_surface_type(&self) -> FixedSurfaceType {
+        self.inner.second_fixed_surface_type()
+    }
+
+    fn second_fixed_surface_value(&self) -> Option<f64> {
+        self.inner.second_fixed_surface_value()
+    }
+
+    fn derived_forecast_type(&self) -> Option<super::tables::DerivedForecastType> {
+        None
+    }
+
+    fn statistical_process_type(&self) -> Option<super::tables::TypeOfStatisticalProcessing> {
+        None
+    }
+
+    fn perturbation_number(&self) -> Option<u8> {
+        Some(self.perturbation_number())
+    }
+
+    fn number_of_ensemble_members(&self) -> Option<u8> {
+        Some(self.number_of_forecasts_in_ensemble())
+    }
+
+    fn wave_period_range(&self) -> Option<WavePeriodRange> {
+        self.inner.wave_period_range()
+    }
+}

gribberish/tests/read.rs

@@ -1697,3 +1697,101 @@ fn read_naqfc_alaska_polar_stereographic_grid() {
     assert_eq!(data.len(), 456225);
     assert!((data[1000] - 31.64).abs() < 0.001, "data[1000]");
 }
+
+/// Summary of the defined (non-missing) values of a decoded field:
+/// `(count, min, max, mean)`.
+fn defined_value_stats(data: &[f64]) -> (usize, f64, f64, f64) {
+    let defined = data.iter().cloned().filter(|v| v.is_finite());
+    let (count, min, max, sum) = defined.fold(
+        (0usize, f64::INFINITY, f64::NEG_INFINITY, 0.0),
+        |(count, min, max, sum), v| (count + 1, min.min(v), max.max(v), sum + v),
+    );
+    (count, min, max, sum / count as f64)
+}
+
+#[test]
+fn read_ecmwf_ifs_wave_period_range_hres() {
+    // ECMWF IFS HRES h1012 (significant wave height, periods 10-12 s), step 24.
+    // IFS encodes its period bands with PDT 4.104 - the ensemble variant of
+    // 4.103 - even for the deterministic run, which is member 0.
+    // Source: ecmwf-open-data 20260921/12z/ifs/0p25/wave/20260921120000-24h-wave-fc.grib2
+    // bytes 0-803074. Expected values from eccodes (grib_get).
+    let grib_data = read_grib_messages("../test-data/ecmwf-ifs-wave-h1012-hres.grib2");
+    let messages = read_messages(grib_data.as_slice()).collect::<Vec<Message>>();
+    assert_eq!(messages.len(), 1);
+
+    let h1012 = &messages[0];
+    assert_eq!(h1012.product_template_id().unwrap(), 104);
+    assert_eq!(h1012.variable_abbrev().unwrap(), "HTSGW");
+    assert_eq!(h1012.unit().unwrap(), "m");
+    assert_eq!(
+        h1012.wave_period_range().unwrap(),
+        Some((Some(10.0), Some(12.0)))
+    );
+    assert_eq!(h1012.perturbation_number().unwrap(), Some(0));
+    assert_eq!(h1012.number_of_ensemble_members().unwrap(), Some(0));
+    assert_eq!(
+        h1012.forecast_date().unwrap(),
+        Utc.with_ymd_and_hms(2026, 9, 22, 12, 0, 0).unwrap()
+    );
+
+    let key = h1012.key().unwrap();
+    assert!(key.contains(":per10-12s"), "key missing period band: {key}");
+    assert!(key.contains(":ens0"), "key missing perturbation: {key}");
+
+    let data = h1012.data().unwrap();
+    assert_eq!(data.len(), 1440 * 721);
+    let (count, min, max, mean) = defined_value_stats(&data);
+    assert_eq!(count, 665628, "defined value count");
+    assert!((min - 0.001245).abs() < 1e-5, "min was {min}");
+    assert!((max - 4.655786).abs() < 1e-4, "max was {max}");
+    assert!((mean - 0.855334).abs() < 1e-4, "mean was {mean}");
+}
+
+#[test]
+fn read_ecmwf_ifs_wave_period_range_ens() {
+    // ECMWF IFS ENS member 1 h1417 (significant wave height, periods 14-17 s),
+    // step 24, PDT 4.104.
+    // Source: ecmwf-open-data 20260921/12z/ifs/0p25/waef/20260921120000-24h-waef-ef.grib2
+    // bytes 102676537-103434707. Expected values from eccodes (grib_get).
+    let grib_data = read_grib_messages("../test-data/ecmwf-ifs-wave-h1417-ens.grib2");
+    let messages = read_messages(grib_data.as_slice()).collect::<Vec<Message>>();
+    assert_eq!(messages.len(), 1);
+
+    let h1417 = &messages[0];
+    assert_eq!(h1417.product_template_id().unwrap(), 104);
+    assert_eq!(h1417.variable_abbrev().unwrap(), "HTSGW");
+    assert_eq!(h1417.unit().unwrap(), "m");
+    assert_eq!(
+        h1417.wave_period_range().unwrap(),
+        Some((Some(14.0), Some(17.0)))
+    );
+    assert_eq!(h1417.perturbation_number().unwrap(), Some(1));
+    assert_eq!(h1417.number_of_ensemble_members().unwrap(), Some(51));
+
+    let key = h1417.key().unwrap();
+    assert!(key.contains(":per14-17s"), "key missing period band: {key}");
+    assert!(key.contains(":ens1"), "key missing perturbation: {key}");
+
+    let data = h1417.data().unwrap();
+    assert_eq!(data.len(), 1440 * 721);
+    let (count, min, max, mean) = defined_value_stats(&data);
+    assert_eq!(count, 665628, "defined value count");
+    assert!((min - 0.001127).abs() < 1e-5, "min was {min}");
+    assert!((max - 6.071440).abs() < 1e-4, "max was {max}");
+    assert!((mean - 0.594546).abs() < 1e-4, "mean was {mean}");
+}
+
+#[test]
+fn probability_template_unit_is_percent() {
+    // PDT 4.5 probabilities carry the thresholded parameter (here PWAT, kg m-2)
+    // but their values are percentages (0-100).
+    let grib_data = read_grib_messages("../test-data/nbm-pwat-prob-above.grib2");
+    let messages = read_messages(grib_data.as_slice()).collect::<Vec<Message>>();
+
+    for message in &messages {
+        assert_eq!(message.product_template_id().unwrap(), 5);
+        assert_eq!(message.variable_abbrev().unwrap(), "PWAT");
+        assert_eq!(message.unit().unwrap(), "%");
+    }
+}