1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
//! Vegetation Phenology Module
//!
//! Extracts phenological metrics from time series data including growing season
//! detection, peak timing, and amplitude calculations.
use crate::error::{Result, TemporalError};
use crate::timeseries::TimeSeriesRaster;
use chrono::Datelike;
use scirs2_core::ndarray::Array3;
use serde::{Deserialize, Serialize};
use tracing::info;
/// Phenology extraction method
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub enum PhenologyMethod {
/// NDVI-based phenology
NDVI,
/// EVI-based phenology
EVI,
/// Threshold-based detection
Threshold,
/// Derivative-based detection
Derivative,
}
/// Phenological metrics
#[derive(Debug, Clone)]
pub struct PhenologyMetrics {
/// Start of growing season (day of year)
pub season_start: Array3<i32>,
/// End of growing season (day of year)
pub season_end: Array3<i32>,
/// Peak vegetation time (day of year)
pub peak_time: Array3<i32>,
/// Maximum vegetation value
pub peak_value: Array3<f64>,
/// Amplitude (peak - base)
pub amplitude: Array3<f64>,
/// Growing season length (days)
pub season_length: Array3<i32>,
}
impl PhenologyMetrics {
/// Create new phenology metrics
#[must_use]
pub fn new(
season_start: Array3<i32>,
season_end: Array3<i32>,
peak_time: Array3<i32>,
peak_value: Array3<f64>,
amplitude: Array3<f64>,
season_length: Array3<i32>,
) -> Self {
Self {
season_start,
season_end,
peak_time,
peak_value,
amplitude,
season_length,
}
}
}
/// Phenology configuration
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct PhenologyConfig {
/// Method to use
pub method: PhenologyMethod,
/// Threshold for season detection (0-1)
pub threshold: f64,
/// Minimum season length (days)
pub min_season_length: i32,
/// Smoothing window size
pub smoothing_window: Option<usize>,
}
impl Default for PhenologyConfig {
fn default() -> Self {
Self {
method: PhenologyMethod::NDVI,
threshold: 0.3,
min_season_length: 60,
smoothing_window: Some(3),
}
}
}
/// Phenology extractor
pub struct PhenologyExtractor;
impl PhenologyExtractor {
/// Extract phenological metrics
///
/// # Errors
/// Returns error if extraction fails
pub fn extract(ts: &TimeSeriesRaster, config: &PhenologyConfig) -> Result<PhenologyMetrics> {
if ts.len() < 10 {
return Err(TemporalError::insufficient_data(
"Need at least 10 observations for phenology",
));
}
let (height, width, n_bands) = ts
.expected_shape()
.ok_or_else(|| TemporalError::insufficient_data("No shape information"))?;
let mut season_start = Array3::from_elem((height, width, n_bands), -1);
let mut season_end = Array3::from_elem((height, width, n_bands), -1);
let mut peak_time = Array3::from_elem((height, width, n_bands), -1);
let mut peak_value = Array3::zeros((height, width, n_bands));
let mut amplitude = Array3::zeros((height, width, n_bands));
let mut season_length = Array3::zeros((height, width, n_bands));
// Collect day of year for each observation
let doys: Vec<i32> = ts
.iter()
.map(|(_, e)| e.metadata.acquisition_date.ordinal() as i32)
.collect();
// Process each pixel
for i in 0..height {
for j in 0..width {
for k in 0..n_bands {
let values = ts.extract_pixel_timeseries(i, j, k)?;
// Apply smoothing if configured
let smoothed = if let Some(window) = config.smoothing_window {
Self::smooth_values(&values, window)
} else {
values.clone()
};
// Find min and max
let min_val = smoothed
.iter()
.copied()
.min_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
.unwrap_or(0.0);
let max_val = smoothed
.iter()
.copied()
.max_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
.unwrap_or(0.0);
let amp = max_val - min_val;
amplitude[[i, j, k]] = amp;
if amp < 0.01 {
// No significant seasonality
continue;
}
// Calculate threshold value
let threshold_val = min_val + config.threshold * amp;
// Find peak
if let Some((peak_idx, &peak)) =
smoothed.iter().enumerate().max_by(|(_, a), (_, b)| {
a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)
})
{
peak_time[[i, j, k]] = doys[peak_idx];
peak_value[[i, j, k]] = peak;
// Find season start (ascending through threshold before peak)
for idx in 0..peak_idx {
if smoothed[idx] >= threshold_val {
season_start[[i, j, k]] = doys[idx];
break;
}
}
// Find season end (descending through threshold after peak)
for idx in (peak_idx + 1)..smoothed.len() {
if smoothed[idx] <= threshold_val {
season_end[[i, j, k]] = doys[idx];
break;
}
}
// Calculate season length
if season_start[[i, j, k]] >= 0 && season_end[[i, j, k]] >= 0 {
let length = season_end[[i, j, k]] - season_start[[i, j, k]];
if length >= config.min_season_length {
season_length[[i, j, k]] = length;
} else {
// Reset if season too short
season_start[[i, j, k]] = -1;
season_end[[i, j, k]] = -1;
}
}
}
}
}
}
info!("Extracted phenology metrics");
Ok(PhenologyMetrics::new(
season_start,
season_end,
peak_time,
peak_value,
amplitude,
season_length,
))
}
/// Simple moving average smoothing
fn smooth_values(values: &[f64], window: usize) -> Vec<f64> {
let mut smoothed = Vec::with_capacity(values.len());
for i in 0..values.len() {
let start = i.saturating_sub(window / 2);
let end = (i + window / 2 + 1).min(values.len());
let window_values = &values[start..end];
let avg = window_values.iter().sum::<f64>() / window_values.len() as f64;
smoothed.push(avg);
}
smoothed
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::timeseries::{TemporalMetadata, TimeSeriesRaster};
use chrono::{DateTime, NaiveDate};
#[test]
fn test_phenology_extraction() {
let mut ts = TimeSeriesRaster::new();
// Simulate growing season (sinusoidal pattern)
for i in 0..365 {
if i % 10 == 0 {
// Sample every 10 days
let dt = DateTime::from_timestamp(1640995200 + i * 86400, 0).expect("valid");
let date =
NaiveDate::from_ymd_opt(2022, 1, 1).expect("valid") + chrono::Duration::days(i);
let metadata = TemporalMetadata::new(dt, date);
// Sinusoidal pattern: peak around day 180
let angle = (i as f64 / 365.0) * 2.0 * std::f64::consts::PI;
let value = 0.3 + 0.4 * angle.sin(); // NDVI-like values 0.3-0.7
let data = Array3::from_elem((5, 5, 1), value);
ts.add_raster(metadata, data).expect("should add");
}
}
let config = PhenologyConfig::default();
let metrics = PhenologyExtractor::extract(&ts, &config).expect("should extract");
// Check that metrics were computed
assert!(metrics.season_start[[2, 2, 0]] > 0);
assert!(metrics.peak_time[[2, 2, 0]] > 0);
assert!(metrics.amplitude[[2, 2, 0]] > 0.0);
}
}