oxigeo_algorithms/raster/
zonal_stats.rs

1//! Zonal statistics computation
2
3use crate::error::{AlgorithmError, Result};
4use oxigeo_core::buffer::RasterBuffer;
5
6/// Statistics for a single zone
7#[derive(Debug, Clone, Copy, Default)]
8pub struct ZonalStatistics {
9    /// Zone ID
10    pub zone_id: i32,
11    /// Count of valid pixels
12    pub count: u64,
13    /// Sum of values
14    pub sum: f64,
15    /// Mean value
16    pub mean: f64,
17    /// Minimum value
18    pub min: f64,
19    /// Maximum value
20    pub max: f64,
21    /// Standard deviation
22    pub std_dev: f64,
23}
24
25/// Computes zonal statistics
26pub fn compute_zonal_stats(
27    values: &RasterBuffer,
28    zones: &RasterBuffer,
29) -> Result<Vec<ZonalStatistics>> {
30    if values.width() != zones.width() || values.height() != zones.height() {
31        return Err(AlgorithmError::InvalidDimensions {
32            message: "Rasters must have same dimensions",
33            actual: values.width() as usize,
34            expected: zones.width() as usize,
35        });
36    }
37
38    let mut stats_map: std::collections::HashMap<i32, ZonalStatistics> =
39        std::collections::HashMap::new();
40
41    // First pass: collect sums and counts
42    for y in 0..values.height() {
43        for x in 0..values.width() {
44            let zone_raw = zones.get_pixel(x, y).map_err(AlgorithmError::Core)?;
45            let value = values.get_pixel(x, y).map_err(AlgorithmError::Core)?;
46
47            // Skip pixels whose zone id is NoData/NaN. Casting a NaN directly
48            // `as i32` saturates to 0 in Rust, which would silently bucket every
49            // out-of-zone pixel into a spurious zone 0 and contaminate its stats.
50            if zones.is_nodata(zone_raw) || !zone_raw.is_finite() {
51                continue;
52            }
53            // Skip NoData/non-finite values so they do not poison sum/mean/stddev.
54            if values.is_nodata(value) || !value.is_finite() {
55                continue;
56            }
57            let zone_id = zone_raw as i32;
58
59            let stats = stats_map.entry(zone_id).or_insert_with(|| ZonalStatistics {
60                zone_id,
61                min: f64::MAX,
62                max: f64::MIN,
63                ..Default::default()
64            });
65
66            stats.count += 1;
67            stats.sum += value;
68            stats.min = stats.min.min(value);
69            stats.max = stats.max.max(value);
70        }
71    }
72
73    // Compute means
74    for stats in stats_map.values_mut() {
75        stats.mean = stats.sum / stats.count as f64;
76    }
77
78    // Second pass: compute standard deviation (same NoData filtering as pass 1
79    // so the variance accumulator matches the mean/count computed above).
80    for y in 0..values.height() {
81        for x in 0..values.width() {
82            let zone_raw = zones.get_pixel(x, y).map_err(AlgorithmError::Core)?;
83            let value = values.get_pixel(x, y).map_err(AlgorithmError::Core)?;
84
85            if zones.is_nodata(zone_raw) || !zone_raw.is_finite() {
86                continue;
87            }
88            if values.is_nodata(value) || !value.is_finite() {
89                continue;
90            }
91            let zone_id = zone_raw as i32;
92
93            if let Some(stats) = stats_map.get_mut(&zone_id) {
94                let diff = value - stats.mean;
95                stats.std_dev += diff * diff;
96            }
97        }
98    }
99
100    // Finalize standard deviations
101    for stats in stats_map.values_mut() {
102        stats.std_dev = (stats.std_dev / stats.count as f64).sqrt();
103    }
104
105    Ok(stats_map.into_values().collect())
106}
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111    use oxigeo_core::types::RasterDataType;
112
113    #[test]
114    fn test_zonal_stats() {
115        let mut values = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
116        let mut zones = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
117
118        // Fill with test data
119        for y in 0..5 {
120            for x in 0..5 {
121                values.set_pixel(x, y, (x + y) as f64).ok();
122                zones.set_pixel(x, y, (x / 2) as f64).ok();
123            }
124        }
125
126        let result = compute_zonal_stats(&values, &zones);
127        assert!(result.is_ok());
128    }
129
130    #[test]
131    fn test_zonal_stats_excludes_nodata_values() {
132        use oxigeo_core::types::NoDataValue;
133        // Value raster: zone 1 everywhere, all values 4.0 except one NaN NoData.
134        let mut values = RasterBuffer::nodata_filled(
135            4,
136            4,
137            RasterDataType::Float64,
138            NoDataValue::Float(f64::NAN),
139        );
140        let mut zones = RasterBuffer::zeros(4, 4, RasterDataType::Float64);
141        for y in 0..4 {
142            for x in 0..4 {
143                values
144                    .set_pixel(x, y, 4.0)
145                    .expect("set value pixel should succeed in test");
146                zones
147                    .set_pixel(x, y, 1.0)
148                    .expect("set zone pixel should succeed in test");
149            }
150        }
151        values
152            .set_pixel(0, 0, f64::NAN)
153            .expect("set NoData value should succeed in test");
154
155        let result =
156            compute_zonal_stats(&values, &zones).expect("zonal stats should succeed in test");
157        assert_eq!(result.len(), 1);
158        let zone = result[0];
159        assert_eq!(zone.zone_id, 1);
160        // 16 cells - 1 NoData = 15 valid cells, all 4.0.
161        assert_eq!(zone.count, 15);
162        assert!(zone.sum.is_finite(), "sum must not be NaN-poisoned");
163        assert!((zone.sum - 60.0).abs() < 1e-9);
164        assert!((zone.mean - 4.0).abs() < 1e-9);
165        assert!((zone.std_dev - 0.0).abs() < 1e-9);
166    }
167
168    #[test]
169    fn test_zonal_stats_skips_nodata_zone() {
170        use oxigeo_core::types::NoDataValue;
171        // Zone raster with a NoData sentinel; those pixels must NOT collapse into
172        // a spurious zone 0.
173        let mut values = RasterBuffer::zeros(3, 3, RasterDataType::Float64);
174        let mut zones = RasterBuffer::nodata_filled(
175            3,
176            3,
177            RasterDataType::Float64,
178            NoDataValue::Float(f64::NAN),
179        );
180        for y in 0..3 {
181            for x in 0..3 {
182                values
183                    .set_pixel(x, y, 10.0)
184                    .expect("set value pixel should succeed in test");
185                zones
186                    .set_pixel(x, y, 5.0)
187                    .expect("set zone pixel should succeed in test");
188            }
189        }
190        // Mark one zone pixel as NoData (NaN).
191        zones
192            .set_pixel(1, 1, f64::NAN)
193            .expect("set NoData zone should succeed in test");
194
195        let result =
196            compute_zonal_stats(&values, &zones).expect("zonal stats should succeed in test");
197        // Only zone 5 should exist -- no spurious zone 0 from the NaN cast.
198        assert_eq!(result.len(), 1);
199        assert_eq!(result[0].zone_id, 5);
200        assert_eq!(result[0].count, 8);
201    }
202}
oxigeo_algorithms/raster/zonal_stats.rs

oxigeo_algorithms/raster/
zonal_stats.rs
