oxigeo_temporal/change/
detection.rs

1//! Change Detection Algorithms
2//!
3//! Implements various change detection algorithms including simple differencing,
4//! statistical tests, and advanced methods like BFAST and LandTrendr.
5
6use crate::error::{Result, TemporalError};
7use crate::stack::RasterStack;
8use crate::timeseries::TimeSeriesRaster;
9use scirs2_core::ndarray::Array3;
10use serde::{Deserialize, Serialize};
11use tracing::info;
12
13/// Change detection method
14#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
15pub enum ChangeDetectionMethod {
16    /// Simple differencing (end - start)
17    SimpleDifference,
18    /// Absolute change
19    AbsoluteChange,
20    /// Relative change (percentage)
21    RelativeChange,
22    /// Z-score based change
23    ZScore,
24    /// Threshold-based change detection
25    Threshold,
26    /// Cumulative sum (CUSUM)
27    CUSUM,
28    /// BFAST (Breaks For Additive Season and Trend)
29    BFAST,
30    /// LandTrendr approximation
31    LandTrendr,
32}
33
34/// Change detection configuration
35#[derive(Debug, Clone, Serialize, Deserialize)]
36pub struct ChangeDetectionConfig {
37    /// Detection method
38    pub method: ChangeDetectionMethod,
39    /// Threshold for threshold-based detection
40    pub threshold: Option<f64>,
41    /// Minimum change magnitude to report
42    pub min_magnitude: Option<f64>,
43    /// NoData value
44    pub nodata: Option<f64>,
45    /// Confidence level for statistical tests
46    pub confidence_level: Option<f64>,
47}
48
49impl Default for ChangeDetectionConfig {
50    fn default() -> Self {
51        Self {
52            method: ChangeDetectionMethod::SimpleDifference,
53            threshold: Some(0.1),
54            min_magnitude: None,
55            nodata: Some(f64::NAN),
56            confidence_level: Some(0.95),
57        }
58    }
59}
60
61/// Change detection result
62#[derive(Debug, Clone)]
63pub struct ChangeDetectionResult {
64    /// Change magnitude
65    pub magnitude: Array3<f64>,
66    /// Change direction (-1, 0, 1)
67    pub direction: Array3<i8>,
68    /// Timestamp of change (if applicable)
69    pub change_time: Option<Array3<i64>>,
70    /// Confidence/significance
71    pub confidence: Option<Array3<f64>>,
72}
73
74impl ChangeDetectionResult {
75    /// Create new change detection result
76    #[must_use]
77    pub fn new(magnitude: Array3<f64>, direction: Array3<i8>) -> Self {
78        Self {
79            magnitude,
80            direction,
81            change_time: None,
82            confidence: None,
83        }
84    }
85
86    /// Add change timestamps
87    #[must_use]
88    pub fn with_change_time(mut self, change_time: Array3<i64>) -> Self {
89        self.change_time = Some(change_time);
90        self
91    }
92
93    /// Add confidence scores
94    #[must_use]
95    pub fn with_confidence(mut self, confidence: Array3<f64>) -> Self {
96        self.confidence = Some(confidence);
97        self
98    }
99}
100
101/// Change detector
102pub struct ChangeDetector;
103
104impl ChangeDetector {
105    /// Detect changes in time series
106    ///
107    /// # Errors
108    /// Returns error if detection fails
109    pub fn detect(
110        ts: &TimeSeriesRaster,
111        config: &ChangeDetectionConfig,
112    ) -> Result<ChangeDetectionResult> {
113        match config.method {
114            ChangeDetectionMethod::SimpleDifference => Self::simple_difference(ts, config),
115            ChangeDetectionMethod::AbsoluteChange => Self::absolute_change(ts, config),
116            ChangeDetectionMethod::RelativeChange => Self::relative_change(ts, config),
117            ChangeDetectionMethod::ZScore => Self::zscore_change(ts, config),
118            ChangeDetectionMethod::Threshold => Self::threshold_change(ts, config),
119            ChangeDetectionMethod::CUSUM => Self::cusum_change(ts, config),
120            ChangeDetectionMethod::BFAST => Self::bfast_change(ts, config),
121            ChangeDetectionMethod::LandTrendr => Self::landtrendr_change(ts, config),
122        }
123    }
124
125    /// Detect changes in raster stack
126    ///
127    /// # Errors
128    /// Returns error if detection fails
129    pub fn detect_stack(
130        stack: &RasterStack,
131        config: &ChangeDetectionConfig,
132    ) -> Result<ChangeDetectionResult> {
133        match config.method {
134            ChangeDetectionMethod::SimpleDifference => Self::simple_difference_stack(stack, config),
135            ChangeDetectionMethod::AbsoluteChange => Self::absolute_change_stack(stack, config),
136            ChangeDetectionMethod::RelativeChange => Self::relative_change_stack(stack, config),
137            _ => Err(TemporalError::change_detection_error(format!(
138                "Method {:?} not supported for stack",
139                config.method
140            ))),
141        }
142    }
143
144    /// Simple difference between first and last observation
145    fn simple_difference(
146        ts: &TimeSeriesRaster,
147        _config: &ChangeDetectionConfig,
148    ) -> Result<ChangeDetectionResult> {
149        if ts.len() < 2 {
150            return Err(TemporalError::insufficient_data(
151                "Need at least 2 observations",
152            ));
153        }
154
155        let first = ts.get_by_index(0)?;
156        let last = ts.get_by_index(ts.len() - 1)?;
157
158        let first_data = first
159            .data
160            .as_ref()
161            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
162        let last_data = last
163            .data
164            .as_ref()
165            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
166
167        let magnitude = last_data - first_data;
168        let direction = Self::compute_direction(&magnitude);
169
170        info!("Completed simple difference change detection");
171        Ok(ChangeDetectionResult::new(magnitude, direction))
172    }
173
174    /// Absolute change
175    fn absolute_change(
176        ts: &TimeSeriesRaster,
177        _config: &ChangeDetectionConfig,
178    ) -> Result<ChangeDetectionResult> {
179        if ts.len() < 2 {
180            return Err(TemporalError::insufficient_data(
181                "Need at least 2 observations",
182            ));
183        }
184
185        let first = ts.get_by_index(0)?;
186        let last = ts.get_by_index(ts.len() - 1)?;
187
188        let first_data = first
189            .data
190            .as_ref()
191            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
192        let last_data = last
193            .data
194            .as_ref()
195            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
196
197        let magnitude = (last_data - first_data).mapv(f64::abs);
198        let direction = Self::compute_direction(&(last_data - first_data));
199
200        info!("Completed absolute change detection");
201        Ok(ChangeDetectionResult::new(magnitude, direction))
202    }
203
204    /// Relative change (percentage)
205    fn relative_change(
206        ts: &TimeSeriesRaster,
207        _config: &ChangeDetectionConfig,
208    ) -> Result<ChangeDetectionResult> {
209        if ts.len() < 2 {
210            return Err(TemporalError::insufficient_data(
211                "Need at least 2 observations",
212            ));
213        }
214
215        let first = ts.get_by_index(0)?;
216        let last = ts.get_by_index(ts.len() - 1)?;
217
218        let first_data = first
219            .data
220            .as_ref()
221            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
222        let last_data = last
223            .data
224            .as_ref()
225            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
226
227        let magnitude =
228            (last_data - first_data) / first_data.mapv(|v| if v != 0.0 { v } else { 1.0 });
229        let direction = Self::compute_direction(&(last_data - first_data));
230
231        info!("Completed relative change detection");
232        Ok(ChangeDetectionResult::new(magnitude, direction))
233    }
234
235    /// Z-score based change detection
236    fn zscore_change(
237        ts: &TimeSeriesRaster,
238        config: &ChangeDetectionConfig,
239    ) -> Result<ChangeDetectionResult> {
240        if ts.len() < 3 {
241            return Err(TemporalError::insufficient_data(
242                "Need at least 3 observations",
243            ));
244        }
245
246        let (height, width, n_bands) = ts
247            .expected_shape()
248            .ok_or_else(|| TemporalError::insufficient_data("No shape information"))?;
249
250        let mut magnitude = Array3::zeros((height, width, n_bands));
251        let mut direction = Array3::zeros((height, width, n_bands));
252
253        let threshold = config.threshold.unwrap_or(2.0);
254
255        for i in 0..height {
256            for j in 0..width {
257                for k in 0..n_bands {
258                    let values = ts.extract_pixel_timeseries(i, j, k)?;
259
260                    let mean = values.iter().sum::<f64>() / values.len() as f64;
261                    let variance = values.iter().map(|v| (v - mean).powi(2)).sum::<f64>()
262                        / values.len() as f64;
263                    let std_dev = variance.sqrt();
264
265                    if std_dev > 0.0 {
266                        let last_value = values[values.len() - 1];
267                        let z_score = ((last_value - mean) / std_dev).abs();
268
269                        magnitude[[i, j, k]] = z_score;
270                        direction[[i, j, k]] = if z_score > threshold {
271                            if last_value > mean { 1 } else { -1 }
272                        } else {
273                            0
274                        };
275                    }
276                }
277            }
278        }
279
280        info!("Completed Z-score change detection");
281        Ok(ChangeDetectionResult::new(magnitude, direction))
282    }
283
284    /// Threshold-based change detection
285    fn threshold_change(
286        ts: &TimeSeriesRaster,
287        config: &ChangeDetectionConfig,
288    ) -> Result<ChangeDetectionResult> {
289        if ts.len() < 2 {
290            return Err(TemporalError::insufficient_data(
291                "Need at least 2 observations",
292            ));
293        }
294
295        let threshold = config.threshold.ok_or_else(|| {
296            TemporalError::invalid_parameter("threshold", "threshold required for this method")
297        })?;
298
299        let first = ts.get_by_index(0)?;
300        let last = ts.get_by_index(ts.len() - 1)?;
301
302        let first_data = first
303            .data
304            .as_ref()
305            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
306        let last_data = last
307            .data
308            .as_ref()
309            .ok_or_else(|| TemporalError::invalid_input("Data not loaded"))?;
310
311        let diff = last_data - first_data;
312        let magnitude = diff.mapv(|v| if v.abs() > threshold { v.abs() } else { 0.0 });
313        let direction = Self::compute_direction(&diff);
314
315        info!("Completed threshold-based change detection");
316        Ok(ChangeDetectionResult::new(magnitude, direction))
317    }
318
319    /// CUSUM change detection
320    fn cusum_change(
321        ts: &TimeSeriesRaster,
322        config: &ChangeDetectionConfig,
323    ) -> Result<ChangeDetectionResult> {
324        if ts.len() < 3 {
325            return Err(TemporalError::insufficient_data(
326                "Need at least 3 observations",
327            ));
328        }
329
330        let (height, width, n_bands) = ts
331            .expected_shape()
332            .ok_or_else(|| TemporalError::insufficient_data("No shape information"))?;
333
334        let mut magnitude = Array3::zeros((height, width, n_bands));
335        let mut direction = Array3::zeros((height, width, n_bands));
336
337        let threshold = config.threshold.unwrap_or(1.0);
338
339        for i in 0..height {
340            for j in 0..width {
341                for k in 0..n_bands {
342                    let values = ts.extract_pixel_timeseries(i, j, k)?;
343
344                    let mean = values.iter().sum::<f64>() / values.len() as f64;
345
346                    // Calculate CUSUM
347                    let mut cusum_pos: f64 = 0.0;
348                    let mut cusum_neg: f64 = 0.0;
349                    let mut max_cusum: f64 = 0.0;
350
351                    for &value in &values {
352                        cusum_pos = (cusum_pos + (value - mean)).max(0.0);
353                        cusum_neg = (cusum_neg - (value - mean)).max(0.0);
354
355                        max_cusum = max_cusum.max(cusum_pos).max(cusum_neg);
356                    }
357
358                    magnitude[[i, j, k]] = max_cusum;
359                    direction[[i, j, k]] = if max_cusum > threshold {
360                        if cusum_pos > cusum_neg { 1 } else { -1 }
361                    } else {
362                        0
363                    };
364                }
365            }
366        }
367
368        info!("Completed CUSUM change detection");
369        Ok(ChangeDetectionResult::new(magnitude, direction))
370    }
371
372    /// BFAST change detection (simplified approximation)
373    fn bfast_change(
374        ts: &TimeSeriesRaster,
375        config: &ChangeDetectionConfig,
376    ) -> Result<ChangeDetectionResult> {
377        // Real BFAST (Verbesselt et al. 2010): harmonic season+trend OLS fit
378        // plus OLS-MOSUM structural-change test; see `crate::change::bfast`.
379        crate::change::bfast::bfast_detect(ts, config)
380    }
381
382    /// LandTrendr change detection (Kennedy et al. 2010 piecewise-linear segmentation).
383    ///
384    /// Delegates per-pixel work to [`crate::change::landtrendr::landtrendr_segment`].
385    /// The reported magnitude is `|last_vertex - first_vertex|` and the direction
386    /// follows the sign of `last_vertex - first_vertex`.
387    fn landtrendr_change(
388        ts: &TimeSeriesRaster,
389        _config: &ChangeDetectionConfig,
390    ) -> Result<ChangeDetectionResult> {
391        info!("LandTrendr piecewise-linear segmentation");
392
393        if ts.len() < 3 {
394            return Err(TemporalError::insufficient_data(
395                "Need at least 3 observations",
396            ));
397        }
398
399        let (height, width, n_bands) = ts
400            .expected_shape()
401            .ok_or_else(|| TemporalError::insufficient_data("No shape information"))?;
402
403        let mut magnitude = Array3::<f64>::zeros((height, width, n_bands));
404        let mut direction = Array3::<i8>::zeros((height, width, n_bands));
405        let options = crate::change::landtrendr::LandTrendrOptions::default();
406
407        for i in 0..height {
408            for j in 0..width {
409                for k in 0..n_bands {
410                    let values = ts.extract_pixel_timeseries(i, j, k)?;
411                    if values.len() < options.min_observations {
412                        continue;
413                    }
414                    if values.iter().any(|v| !v.is_finite()) {
415                        continue;
416                    }
417                    if let Ok(result) =
418                        crate::change::landtrendr::landtrendr_segment(&values, &options)
419                        && !result.vertices.is_empty()
420                    {
421                        let first = result.vertices.first().map(|v| v.value).unwrap_or(0.0);
422                        let last = result.vertices.last().map(|v| v.value).unwrap_or(0.0);
423                        magnitude[[i, j, k]] = (last - first).abs();
424                        direction[[i, j, k]] = if last > first {
425                            1
426                        } else if last < first {
427                            -1
428                        } else {
429                            0
430                        };
431                    }
432                }
433            }
434        }
435
436        info!("Completed LandTrendr change detection");
437        Ok(ChangeDetectionResult::new(magnitude, direction))
438    }
439
440    /// Simple difference for raster stack
441    fn simple_difference_stack(
442        stack: &RasterStack,
443        _config: &ChangeDetectionConfig,
444    ) -> Result<ChangeDetectionResult> {
445        let (n_time, _height, _width, _n_bands) = stack.shape();
446
447        if n_time < 2 {
448            return Err(TemporalError::insufficient_data(
449                "Need at least 2 time steps",
450            ));
451        }
452
453        let data = stack.data();
454        let first = data.slice(s![0, .., .., ..]).to_owned();
455        let last = data.slice(s![n_time - 1, .., .., ..]).to_owned();
456
457        let magnitude = &last - &first;
458        let direction = Self::compute_direction(&magnitude);
459
460        Ok(ChangeDetectionResult::new(magnitude, direction))
461    }
462
463    /// Absolute change for raster stack
464    fn absolute_change_stack(
465        stack: &RasterStack,
466        _config: &ChangeDetectionConfig,
467    ) -> Result<ChangeDetectionResult> {
468        let (n_time, _height, _width, _n_bands) = stack.shape();
469
470        if n_time < 2 {
471            return Err(TemporalError::insufficient_data(
472                "Need at least 2 time steps",
473            ));
474        }
475
476        let data = stack.data();
477        let first = data.slice(s![0, .., .., ..]).to_owned();
478        let last = data.slice(s![n_time - 1, .., .., ..]).to_owned();
479
480        let diff = &last - &first;
481        let magnitude = diff.mapv(f64::abs);
482        let direction = Self::compute_direction(&diff);
483
484        Ok(ChangeDetectionResult::new(magnitude, direction))
485    }
486
487    /// Relative change for raster stack
488    fn relative_change_stack(
489        stack: &RasterStack,
490        _config: &ChangeDetectionConfig,
491    ) -> Result<ChangeDetectionResult> {
492        let (n_time, _height, _width, _n_bands) = stack.shape();
493
494        if n_time < 2 {
495            return Err(TemporalError::insufficient_data(
496                "Need at least 2 time steps",
497            ));
498        }
499
500        let data = stack.data();
501        let first = data.slice(s![0, .., .., ..]).to_owned();
502        let last = data.slice(s![n_time - 1, .., .., ..]).to_owned();
503
504        let diff = &last - &first;
505        let magnitude = &diff / &first.mapv(|v| if v != 0.0 { v } else { 1.0 });
506        let direction = Self::compute_direction(&diff);
507
508        Ok(ChangeDetectionResult::new(magnitude, direction))
509    }
510
511    /// Compute change direction
512    fn compute_direction(change: &Array3<f64>) -> Array3<i8> {
513        let shape = change.shape();
514        let mut direction = Array3::zeros((shape[0], shape[1], shape[2]));
515
516        for i in 0..shape[0] {
517            for j in 0..shape[1] {
518                for k in 0..shape[2] {
519                    let c = change[[i, j, k]];
520                    direction[[i, j, k]] = if c > 0.0 {
521                        1
522                    } else if c < 0.0 {
523                        -1
524                    } else {
525                        0
526                    };
527                }
528            }
529        }
530
531        direction
532    }
533}
534
535// Import ndarray slice macro
536use scirs2_core::ndarray::s;
537
538#[cfg(test)]
539mod tests {
540    use super::*;
541    use crate::timeseries::{TemporalMetadata, TimeSeriesRaster};
542    use chrono::{DateTime, NaiveDate};
543
544    #[test]
545    fn test_simple_difference() {
546        let mut ts = TimeSeriesRaster::new();
547
548        let dt1 = DateTime::from_timestamp(1640995200, 0).expect("valid");
549        let date1 = NaiveDate::from_ymd_opt(2022, 1, 1).expect("valid");
550        let metadata1 = TemporalMetadata::new(dt1, date1);
551        ts.add_raster(metadata1, Array3::from_elem((5, 5, 1), 10.0))
552            .expect("should add");
553
554        let dt2 = DateTime::from_timestamp(1641081600, 0).expect("valid");
555        let date2 = NaiveDate::from_ymd_opt(2022, 1, 2).expect("valid");
556        let metadata2 = TemporalMetadata::new(dt2, date2);
557        ts.add_raster(metadata2, Array3::from_elem((5, 5, 1), 20.0))
558            .expect("should add");
559
560        let config = ChangeDetectionConfig::default();
561        let result = ChangeDetector::detect(&ts, &config).expect("should detect");
562
563        assert_eq!(result.magnitude[[0, 0, 0]], 10.0);
564        assert_eq!(result.direction[[0, 0, 0]], 1);
565    }
566}
oxigeo_temporal/change/detection.rs

oxigeo_temporal/change/
detection.rs
