From a9eed7f3784ec8d22574481e6a8cb33a8c429be1 Mon Sep 17 00:00:00 2001 From: Schuwi Date: Fri, 6 Mar 2026 21:40:39 +0100 Subject: [PATCH] Fix GRIB2 scale factor decoding (sign-magnitude, not two's complement) GRIB2 encodes binary (E) and decimal (D) scale factors in sign-and-magnitude format (MSB = sign bit, lower 15 bits = magnitude), not two's complement. The previous read_i16 call misinterpreted 0x8009 as -32759 instead of -9, causing 2^E to underflow to 0.0 in f32 and zeroing all decoded values. Fix: replace read_i16 with read_grib_scale for E and D in decode_simple_packing. Add regression tests covering sign-magnitude decoding and a realistic DWD CLCT packing scenario. --- src/grib.rs | 52 +++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/src/grib.rs b/src/grib.rs index 58ec429..fe64f05 100644 --- a/src/grib.rs +++ b/src/grib.rs @@ -171,20 +171,23 @@ fn decode_simple_packing(sec5: &[u8], sec7: &[u8], n: u32) -> Result> { let r_bytes: [u8; 4] = sec5[0..4].try_into().unwrap(); let r = f32::from_bits(u32::from_be_bytes(r_bytes)); - let e = read_i16(&sec5[4..6]); - let d = read_i16(&sec5[6..8]); + // GRIB2 encodes E and D as sign-and-magnitude (MSB = sign, lower 15 bits = magnitude), + // NOT two's complement. Using i16::from_be_bytes would misinterpret e.g. 0x8009 as + // -32759 instead of the correct -9, causing 2^E to underflow to zero. + let e = read_grib_scale(&sec5[4..6]); + let d = read_grib_scale(&sec5[6..8]); let n_bits = sec5[8] as u32; if n_bits == 0 { // All values equal the reference value - return Ok(vec![r * 10f32.powi(-(d as i32)); n as usize]); + return Ok(vec![r * 10f32.powi(-d); n as usize]); } if n_bits > 32 { bail!("Implausible n_bits = {}", n_bits); } - let scale_e = 2f32.powi(e as i32); - let scale_d = 10f32.powi(-(d as i32)); + let scale_e = 2f32.powi(e); + let scale_d = 10f32.powi(-d); let mut values = Vec::with_capacity(n as usize); let mut bit_offset: u64 = 0; @@ -239,8 +242,15 @@ fn read_u16(b: &[u8]) -> u16 { u16::from_be_bytes(b[0..2].try_into().unwrap()) } -fn read_i16(b: &[u8]) -> i16 { - i16::from_be_bytes(b[0..2].try_into().unwrap()) +/// Decode a GRIB2 sign-and-magnitude 16-bit integer. +/// +/// The WMO GRIB2 specification encodes binary and decimal scale factors using +/// a non-standard representation: bit 15 is the sign (1 = negative) and bits +/// 14–0 hold the unsigned magnitude. This is different from two's complement. +fn read_grib_scale(b: &[u8]) -> i32 { + let raw = u16::from_be_bytes(b[0..2].try_into().unwrap()); + let magnitude = (raw & 0x7FFF) as i32; + if raw & 0x8000 != 0 { -magnitude } else { magnitude } } #[cfg(test)] @@ -270,6 +280,34 @@ mod tests { assert_eq!(result, vec![5.0f32, 5.0, 5.0]); } + #[test] + fn test_read_grib_scale_sign_magnitude() { + // 0x8009 = sign bit set, magnitude 9 → should be -9, not -32759 (two's complement) + assert_eq!(read_grib_scale(&[0x80, 0x09]), -9); + // 0x0009 = positive 9 + assert_eq!(read_grib_scale(&[0x00, 0x09]), 9); + // 0x0000 = zero + assert_eq!(read_grib_scale(&[0x00, 0x00]), 0); + // 0x8000 = negative zero (sign bit set, magnitude 0) → 0 + assert_eq!(read_grib_scale(&[0x80, 0x00]), 0); + } + + #[test] + fn test_simple_packing_with_neg_e() { + // Simulate DWD CLCT packing: R=0, E=-9 (0x8009), D=0, n_bits=16 + // A packed value of 51200 should decode to 51200 * 2^-9 = 100.0 + let r = 0.0f32; + let mut sec5 = vec![0u8; 10]; + sec5[0..4].copy_from_slice(&r.to_bits().to_be_bytes()); + sec5[4..6].copy_from_slice(&[0x80, 0x09]); // E = -9 in sign-magnitude + sec5[6..8].copy_from_slice(&[0x00, 0x00]); // D = 0 + sec5[8] = 16; // 16 bits per value + // Packed value 51200 as big-endian u16 + let sec7 = [0xC8u8, 0x00]; // 0xC800 = 51200 + let result = decode_simple_packing(&sec5, &sec7, 1).unwrap(); + assert!((result[0] - 100.0f32).abs() < 0.01, "Expected ~100.0, got {}", result[0]); + } + #[test] fn test_expand_with_bitmap() { // bitmap: 0b10110000 → positions 0, 2, 3 are present; 1, 4-7 are missing