oxigdal_algorithms/raster/
zonal_stats.rs

1//! Zonal statistics computation
2
3use crate::error::{AlgorithmError, Result};
4use oxigdal_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_id = zones.get_pixel(x, y).map_err(AlgorithmError::Core)? as i32;
45            let value = values.get_pixel(x, y).map_err(AlgorithmError::Core)?;
46
47            let stats = stats_map.entry(zone_id).or_insert_with(|| ZonalStatistics {
48                zone_id,
49                min: f64::MAX,
50                max: f64::MIN,
51                ..Default::default()
52            });
53
54            stats.count += 1;
55            stats.sum += value;
56            stats.min = stats.min.min(value);
57            stats.max = stats.max.max(value);
58        }
59    }
60
61    // Compute means
62    for stats in stats_map.values_mut() {
63        stats.mean = stats.sum / stats.count as f64;
64    }
65
66    // Second pass: compute standard deviation
67    for y in 0..values.height() {
68        for x in 0..values.width() {
69            let zone_id = zones.get_pixel(x, y).map_err(AlgorithmError::Core)? as i32;
70            let value = values.get_pixel(x, y).map_err(AlgorithmError::Core)?;
71
72            if let Some(stats) = stats_map.get_mut(&zone_id) {
73                let diff = value - stats.mean;
74                stats.std_dev += diff * diff;
75            }
76        }
77    }
78
79    // Finalize standard deviations
80    for stats in stats_map.values_mut() {
81        stats.std_dev = (stats.std_dev / stats.count as f64).sqrt();
82    }
83
84    Ok(stats_map.into_values().collect())
85}
86
87#[cfg(test)]
88mod tests {
89    use super::*;
90    use oxigdal_core::types::RasterDataType;
91
92    #[test]
93    fn test_zonal_stats() {
94        let mut values = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
95        let mut zones = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
96
97        // Fill with test data
98        for y in 0..5 {
99            for x in 0..5 {
100                values.set_pixel(x, y, (x + y) as f64).ok();
101                zones.set_pixel(x, y, (x / 2) as f64).ok();
102            }
103        }
104
105        let result = compute_zonal_stats(&values, &zones);
106        assert!(result.is_ok());
107    }
108}
oxigdal_algorithms/raster/zonal_stats.rs

oxigdal_algorithms/raster/
zonal_stats.rs
