Add PDT 107 support for anomaly products with reference to normal (#139)

PDT 107 (Derived Forecasts with Reference to a Normal) was silently
dropped, causing 16 anomaly messages in S2S statistical files to be
missing. This adds a new template with the correct byte layout (PDT 8
time-interval offsets + derived forecast at byte 58), an `is_anomaly`
trait method, and `:anom` key suffix to prevent collisions with PDT 12
messages sharing the same variable/level/stat/derived type.

Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: jbeezley

gribberish/src/message.rs

@@ -229,8 +229,14 @@ impl<'a> Message<'a> {
             None => "".to_string(),
         };
 
+        let anomaly = if self.is_anomaly().unwrap_or(false) {
+            ":anom"
+        } else {
+            ""
+        };
+
         Ok(format!(
-            "{var}{time}{first_level}{second_level}{perturbation}{percentile}{probability}:{statistical_process}{generating_process}{derived_forecast_type}"
+            "{var}{time}{first_level}{second_level}{perturbation}{percentile}{probability}{anomaly}:{statistical_process}{generating_process}{derived_forecast_type}"
         ))
     }
 
@@ -532,6 +538,16 @@ impl<'a> Message<'a> {
         }
     }
 
+    pub fn is_anomaly(&self) -> Result<bool, GribberishError> {
+        match self {
+            Message::Grib1 { .. } => Ok(false),
+            Message::Grib2 { .. } => {
+                let product_template = self.product_template()?;
+                Ok(product_template.is_anomaly())
+            }
+        }
+    }
+
     pub fn perturbation_number(&self) -> Result<Option<u8>, GribberishError> {
         match self {
             Message::Grib1 { .. } => Ok(None),

gribberish/src/message_metadata.rs

@@ -53,6 +53,7 @@ pub struct MessageMetadata {
     pub forecast_probability_number: Option<u8>,
     pub probability_lower_limit: Option<f64>,
     pub probability_upper_limit: Option<f64>,
+    pub is_anomaly: bool,
 }
 
 impl MessageMetadata {
@@ -175,6 +176,7 @@ impl<'a> TryFrom<&Message<'a>> for MessageMetadata {
             forecast_probability_number: message.forecast_probability_number()?,
             probability_lower_limit: message.probability_lower_limit()?,
             probability_upper_limit: message.probability_upper_limit()?,
+            is_anomaly: message.is_anomaly()?,
         })
     }
 }

gribberish/src/parser.rs

@@ -329,6 +329,7 @@ fn grib1_to_metadata(message: &Grib1Message) -> Result<MessageMetadata, Gribberi
         forecast_probability_number: None,
         probability_lower_limit: None,
         probability_upper_limit: None,
+        is_anomaly: false,
     })
 }
 

gribberish/src/sections/product_definition.rs

@@ -1,6 +1,7 @@
 use super::grib_section::GribSection;
 use crate::{
     templates::product::{
+        derived_ensemble_forecast_time_interval_reference_template::DerivedEnsembleForecastTimeIntervalReferenceTemplate,
         derived_ensemble_horizontal_forecast_time_interval_template::DerivedEnsembleHorizontalForecastTimeIntervalTemplate,
         product_template::ProductTemplate,
         AverageAccumulationExtremeHorizontalAnalysisForecastTemplate,
@@ -73,6 +74,12 @@ impl<'a> ProductDefinitionSection<'a> {
                     discipline,
                 ),
             )),
+            107 => Some(Box::new(
+                DerivedEnsembleForecastTimeIntervalReferenceTemplate::new(
+                    self.data.to_vec(),
+                    discipline,
+                ),
+            )),
             _ => None,
         }
     }

gribberish/src/templates/product/derived_ensemble_forecast_time_interval_reference_template.rs

@@ -0,0 +1,171 @@
+use crate::templates::template::{Template, TemplateType};
+use crate::utils::{read_u16_from_bytes, read_u32_from_bytes};
+use chrono::{prelude::*, Duration};
+
+use super::product_template::ProductTemplate;
+use super::tables::{
+    DerivedForecastType, FixedSurfaceType, GeneratingProcess, TimeUnit,
+    TypeOfStatisticalProcessing, TypeOfTimeInterval,
+};
+use super::HorizontalAnalysisForecastTemplate;
+
+pub struct DerivedEnsembleForecastTimeIntervalReferenceTemplate {
+    data: Vec<u8>,
+    discipline: u8,
+}
+
+impl Template for DerivedEnsembleForecastTimeIntervalReferenceTemplate {
+    fn data(&self) -> &[u8] {
+        &self.data
+    }
+
+    fn template_number(&self) -> u16 {
+        107
+    }
+
+    fn template_type(&self) -> TemplateType {
+        TemplateType::Product
+    }
+
+    fn template_name(&self) -> &str {
+        "Derived forecasts based on all ensemble members at a horizontal level
+        or in a horizontal layer in a continuous or non-continuous time interval
+        with reference to a normal"
+    }
+}
+
+impl DerivedEnsembleForecastTimeIntervalReferenceTemplate {
+    pub fn new(data: Vec<u8>, discipline: u8) -> Self {
+        Self { data, discipline }
+    }
+
+    pub fn first_fixed_surface_scale_factor(&self) -> i8 {
+        as_signed!(self.data[23], 8, i8)
+    }
+
+    pub fn first_fixed_surface_scaled_value(&self) -> i32 {
+        as_signed!(read_u32_from_bytes(&self.data, 24).unwrap_or(0), 32, i32)
+    }
+
+    pub fn second_fixed_surface_scale_factor(&self) -> i8 {
+        as_signed!(self.data[29], 8, i8)
+    }
+
+    pub fn second_fixed_surface_scaled_value(&self) -> i32 {
+        as_signed!(read_u32_from_bytes(&self.data, 30).unwrap_or(0), 32, i32)
+    }
+
+    pub fn valid_end_date(&self) -> DateTime<Utc> {
+        let data = self.data();
+        let year = read_u16_from_bytes(data, 34).unwrap_or(0) as i32;
+        let month = data[36] as u32;
+        let day = data[37] as u32;
+        let hour = data[38] as u32;
+        let minute = data[39] as u32;
+        let second = data[40] as u32;
+
+        Utc.with_ymd_and_hms(year as i32, month, day, hour, minute, second)
+            .unwrap()
+    }
+
+    pub fn number_of_time_ranges(&self) -> u8 {
+        self.data()[41]
+    }
+
+    pub fn number_of_values_missing_from_stats(&self) -> u32 {
+        read_u32_from_bytes(self.data(), 42).unwrap_or(0)
+    }
+
+    pub fn type_of_time_interval(&self) -> TypeOfTimeInterval {
+        self.data()[47].into()
+    }
+
+    pub fn statistical_process_time_unit(&self) -> TimeUnit {
+        self.data()[48].into()
+    }
+
+    pub fn statistical_process_time_interval(&self) -> u32 {
+        read_u32_from_bytes(self.data(), 49).unwrap_or(0)
+    }
+
+    pub fn number_of_forecasts_in_ensemble(&self) -> u8 {
+        self.data[59]
+    }
+}
+
+impl ProductTemplate for DerivedEnsembleForecastTimeIntervalReferenceTemplate {
+    fn discipline(&self) -> u8 {
+        self.discipline
+    }
+
+    fn category_value(&self) -> u8 {
+        self.data[9]
+    }
+
+    fn parameter_value(&self) -> u8 {
+        self.data[10]
+    }
+
+    fn generating_process(&self) -> GeneratingProcess {
+        self.data[11].into()
+    }
+
+    fn time_unit(&self) -> TimeUnit {
+        self.data[17].into()
+    }
+
+    fn time_increment_unit(&self) -> Option<TimeUnit> {
+        Some(self.data()[53].into())
+    }
+
+    fn time_interval(&self) -> u32 {
+        read_u32_from_bytes(&self.data, 18).unwrap_or(0)
+    }
+
+    fn time_increment_interval(&self) -> Option<u32> {
+        Some(read_u32_from_bytes(self.data(), 54).unwrap_or(0))
+    }
+
+    fn forecast_datetime(&self, reference_date: DateTime<Utc>) -> DateTime<Utc> {
+        let offset_duration: Duration = self.time_interval_duration();
+        reference_date + offset_duration
+    }
+
+    fn forecast_end_datetime(&self, _reference_date: DateTime<Utc>) -> Option<DateTime<Utc>> {
+        Some(self.valid_end_date())
+    }
+
+    fn first_fixed_surface_type(&self) -> FixedSurfaceType {
+        self.data[22].into()
+    }
+
+    fn first_fixed_surface_value(&self) -> Option<f64> {
+        HorizontalAnalysisForecastTemplate::scale_value(
+            self.first_fixed_surface_scale_factor(),
+            self.first_fixed_surface_scaled_value(),
+        )
+    }
+
+    fn second_fixed_surface_type(&self) -> FixedSurfaceType {
+        self.data[28].into()
+    }
+
+    fn second_fixed_surface_value(&self) -> Option<f64> {
+        HorizontalAnalysisForecastTemplate::scale_value(
+            self.second_fixed_surface_scale_factor(),
+            self.second_fixed_surface_scaled_value(),
+        )
+    }
+
+    fn derived_forecast_type(&self) -> Option<DerivedForecastType> {
+        Some(self.data[58].into())
+    }
+
+    fn statistical_process_type(&self) -> Option<TypeOfStatisticalProcessing> {
+        Some(self.data()[46].into())
+    }
+
+    fn is_anomaly(&self) -> bool {
+        true
+    }
+}

gribberish/src/templates/product/mod.rs

@@ -1,4 +1,5 @@
 pub mod average_accumulation_extreme_horizontal_analysis_template;
+pub mod derived_ensemble_forecast_time_interval_reference_template;
 pub mod derived_ensemble_horizontal_analysis_template;
 pub mod derived_ensemble_horizontal_forecast_time_interval_template;
 pub mod ensemble_forecast_time_interval_template;
@@ -12,6 +13,7 @@ pub mod product_template;
 pub mod tables;
 
 pub use average_accumulation_extreme_horizontal_analysis_template::AverageAccumulationExtremeHorizontalAnalysisForecastTemplate;
+pub use derived_ensemble_forecast_time_interval_reference_template::DerivedEnsembleForecastTimeIntervalReferenceTemplate;
 pub use derived_ensemble_horizontal_analysis_template::DerivedEnsembleHorizontalAnalysisForecastTemplate;
 pub use ensemble_forecast_time_interval_template::EnsembleForecastTimeIntervalTemplate;
 pub use horizontal_analysis_template::HorizontalAnalysisForecastTemplate;

gribberish/src/templates/product/product_template.rs

@@ -59,6 +59,10 @@ pub trait ProductTemplate {
         None
     }
 
+    fn is_anomaly(&self) -> bool {
+        false
+    }
+
     fn category(&self) -> &'static str {
         category(self.discipline(), self.category_value())
     }

gribberish/tests/read.rs

@@ -1,6 +1,7 @@
 extern crate gribberish;
 
 use gribberish::message::{read_messages, Message};
+use gribberish::templates::product::tables::{DerivedForecastType, TypeOfStatisticalProcessing};
 use std::time::Instant;
 use std::vec::Vec;
 
@@ -977,3 +978,89 @@ fn read_percentile_and_probability_templates() {
         dups
     );
 }
+
+#[test]
+fn read_pdt107_anomaly_with_reference() {
+    // Fixture: 4 messages with zeroed data, 0.5° global grid.
+    // Messages 0-1: PDT 12 (derived ensemble time interval)
+    // Messages 2-3: PDT 107 (same, with reference to normal — anomaly)
+    // Messages 0 and 2 share the same variable/level/stat/derived type
+    // and would collide without the :anom key suffix.
+    let grib_data = read_grib_messages("../test-data/s2s-pdt12-pdt107-anomaly.grib2");
+    let messages = read_messages(grib_data.as_slice()).collect::<Vec<Message>>();
+
+    assert_eq!(messages.len(), 4, "Expected 4 messages in fixture");
+
+    // --- PDT 12: TMP 2m Max UnweightedMean ---
+    let pdt12_max = &messages[0];
+    assert_eq!(pdt12_max.product_template_id().unwrap(), 12);
+    assert_eq!(pdt12_max.variable_abbrev().unwrap(), "TMP");
+    assert!(!pdt12_max.is_anomaly().unwrap());
+
+    // --- PDT 12: TMP 2m Avg UnweightedMean ---
+    let pdt12_avg = &messages[1];
+    assert_eq!(pdt12_avg.product_template_id().unwrap(), 12);
+    assert_eq!(pdt12_avg.variable_abbrev().unwrap(), "TMP");
+    assert!(!pdt12_avg.is_anomaly().unwrap());
+
+    // --- PDT 107: TMP 2m Max UnweightedMean (anomaly) ---
+    let pdt107_max = &messages[2];
+    assert_eq!(pdt107_max.product_template_id().unwrap(), 107);
+    assert_eq!(pdt107_max.variable_abbrev().unwrap(), "TMP");
+    assert!(pdt107_max.is_anomaly().unwrap());
+    assert_eq!(
+        pdt107_max.derived_forecast_type().unwrap(),
+        Some(DerivedForecastType::UnweightedMean)
+    );
+    assert_eq!(
+        pdt107_max.statistical_process_type().unwrap(),
+        Some(TypeOfStatisticalProcessing::Maximum)
+    );
+    let data = pdt107_max.data().unwrap();
+    assert_eq!(data.len(), 259920);
+
+    // --- PDT 107: TMP 1000hPa Avg UnweightedMean (anomaly) ---
+    let pdt107_avg = &messages[3];
+    assert_eq!(pdt107_avg.product_template_id().unwrap(), 107);
+    assert_eq!(pdt107_avg.variable_abbrev().unwrap(), "TMP");
+    assert!(pdt107_avg.is_anomaly().unwrap());
+    assert_eq!(
+        pdt107_avg.derived_forecast_type().unwrap(),
+        Some(DerivedForecastType::UnweightedMean)
+    );
+    assert_eq!(
+        pdt107_avg.statistical_process_type().unwrap(),
+        Some(TypeOfStatisticalProcessing::Average)
+    );
+
+    // --- All 4 messages should have unique keys ---
+    let mut keys = Vec::new();
+    let mut dups = Vec::new();
+    for message in &messages {
+        let key = message.key().unwrap();
+        if keys.contains(&key) {
+            dups.push(key);
+        } else {
+            keys.push(key);
+        }
+    }
+    assert_eq!(
+        dups.len(),
+        0,
+        "Found {} duplicate keys: {:?}",
+        dups.len(),
+        dups
+    );
+
+    // Verify the anomaly key contains :anom and the non-anomaly does not
+    let pdt12_key = pdt12_max.key().unwrap();
+    let pdt107_key = pdt107_max.key().unwrap();
+    assert!(
+        !pdt12_key.contains(":anom"),
+        "PDT 12 key should not contain :anom"
+    );
+    assert!(
+        pdt107_key.contains(":anom"),
+        "PDT 107 key should contain :anom"
+    );
+}

python/src/dataset.rs

@@ -156,8 +156,9 @@ pub fn parse_grib_dataset<'py>(
             }
             _ => "".to_string(),
         };
+        let anom = if v.2.is_anomaly { "_anom" } else { "" };
         let hash = format!(
-            "{surf}_{stat}{gen}{accum_period}{derived}{prob_type}{prob_limits}",
+            "{surf}_{stat}{gen}{accum_period}{derived}{prob_type}{prob_limits}{anom}",
             surf = v.2.first_fixed_surface_type.coordinate_name(),
             stat =
                 v.2.statistical_process