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

@@ -171,20 +171,23 @@ fn decode_simple_packing(sec5: &[u8], sec7: &[u8], n: u32) -> Result<Vec<f32>> {
 
     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