use anyhow::{Context, Result};
use std::path::PathBuf;
use std::time::Instant;
use surtgis_algorithms::imagery::{RasterDiffParams, raster_difference};
use surtgis_algorithms::temporal::{
AnomalyMethod, PhenologyParams, linear_trend, mann_kendall, sens_slope, temporal_anomaly,
temporal_percentile, temporal_stats, vegetation_phenology,
};
use crate::commands::TemporalCommands;
use crate::helpers::{done, read_dem, write_result};
pub fn handle(cmd: TemporalCommands, compress: bool) -> Result<()> {
match cmd {
TemporalCommands::Stats {
input,
outdir,
stats,
} => handle_stats(input, outdir, stats, compress),
TemporalCommands::Trend {
input,
outdir,
method,
times,
} => handle_trend(input, outdir, method, times, compress),
TemporalCommands::Change {
before,
after,
output,
decrease_threshold,
increase_threshold,
} => handle_change(
before,
after,
output,
decrease_threshold,
increase_threshold,
compress,
),
TemporalCommands::Anomaly {
reference,
target,
outdir,
method,
} => handle_anomaly(reference, target, outdir, method, compress),
TemporalCommands::Phenology {
input,
outdir,
doys,
threshold,
smooth,
} => handle_phenology(input, outdir, doys, threshold, smooth, compress),
}
}
fn load_raster_stack(paths: &[PathBuf]) -> Result<Vec<surtgis_core::Raster<f64>>> {
println!("Loading {} rasters...", paths.len());
let mut rasters = Vec::with_capacity(paths.len());
for (i, path) in paths.iter().enumerate() {
let r = read_dem(path)
.with_context(|| format!("Failed to read raster {}: {}", i, path.display()))?;
rasters.push(r);
}
Ok(rasters)
}
fn parse_times(times_str: &str) -> Result<Vec<f64>> {
times_str
.split(',')
.map(|s| {
s.trim()
.parse::<f64>()
.with_context(|| format!("invalid time value: '{}'", s))
})
.collect()
}
fn handle_stats(
input: Vec<PathBuf>,
outdir: PathBuf,
stats_str: String,
compress: bool,
) -> Result<()> {
let rasters = load_raster_stack(&input)?;
let refs: Vec<&surtgis_core::Raster<f64>> = rasters.iter().collect();
std::fs::create_dir_all(&outdir)?;
let start = Instant::now();
let requested: Vec<&str> = stats_str.split(',').map(|s| s.trim()).collect();
let has_basic = requested
.iter()
.any(|s| matches!(*s, "mean" | "std" | "min" | "max" | "count"));
if has_basic {
println!("Computing temporal statistics...");
let ts = temporal_stats(&refs).context("temporal_stats failed")?;
for stat_name in &requested {
let (raster, name) = match *stat_name {
"mean" => (Some(&ts.mean), "mean"),
"std" => (Some(&ts.std), "std"),
"min" => (Some(&ts.min), "min"),
"max" => (Some(&ts.max), "max"),
"count" => (Some(&ts.count), "count"),
_ => (None, *stat_name),
};
if let Some(r) = raster {
let path = outdir.join(format!("{}.tif", name));
write_result(r, &path, compress)?;
println!(" {} → {}", name, path.display());
}
}
}
for stat_name in &requested {
if stat_name.starts_with('p') {
if let Ok(pct) = stat_name[1..].parse::<f64>() {
println!("Computing percentile {}...", pct);
let r = temporal_percentile(&refs, pct)
.with_context(|| format!("percentile {} failed", pct))?;
let path = outdir.join(format!("p{}.tif", pct as u32));
write_result(&r, &path, compress)?;
println!(" p{} → {}", pct as u32, path.display());
}
}
}
done("Temporal statistics", &outdir, start.elapsed());
Ok(())
}
fn handle_trend(
input: Vec<PathBuf>,
outdir: PathBuf,
method: String,
times_str: Option<String>,
compress: bool,
) -> Result<()> {
let rasters = load_raster_stack(&input)?;
let refs: Vec<&surtgis_core::Raster<f64>> = rasters.iter().collect();
let times = times_str.as_ref().map(|s| parse_times(s)).transpose()?;
std::fs::create_dir_all(&outdir)?;
let start = Instant::now();
match method.as_str() {
"linear" | "ols" => {
println!("Computing linear trend (OLS)...");
let result = linear_trend(&refs, times.as_deref()).context("linear_trend failed")?;
let slope_path = outdir.join("slope.tif");
let intercept_path = outdir.join("intercept.tif");
let r2_path = outdir.join("r_squared.tif");
let pval_path = outdir.join("p_value.tif");
write_result(&result.slope, &slope_path, compress)?;
write_result(&result.intercept, &intercept_path, compress)?;
write_result(&result.r_squared, &r2_path, compress)?;
write_result(&result.p_value, &pval_path, compress)?;
println!(" slope → {}", slope_path.display());
println!(" intercept → {}", intercept_path.display());
println!(" R² → {}", r2_path.display());
println!(" p-value → {}", pval_path.display());
}
"mann-kendall" | "mk" => {
println!("Computing Mann-Kendall trend test...");
let result = mann_kendall(&refs).context("mann_kendall failed")?;
let tau_path = outdir.join("tau.tif");
let pval_path = outdir.join("p_value.tif");
let trend_path = outdir.join("trend.tif");
let sens_path = outdir.join("sens_slope.tif");
write_result(&result.tau, &tau_path, compress)?;
write_result(&result.p_value, &pval_path, compress)?;
write_result(&result.trend, &trend_path, compress)?;
write_result(&result.sens_slope, &sens_path, compress)?;
println!(" tau → {}", tau_path.display());
println!(" p-value → {}", pval_path.display());
println!(
" trend → {} (1=up, 0=none, -1=down)",
trend_path.display()
);
println!(" Sen's slope → {}", sens_path.display());
}
"sens" | "sens-slope" => {
println!("Computing Sen's slope...");
let result = sens_slope(&refs).context("sens_slope failed")?;
let path = outdir.join("sens_slope.tif");
write_result(&result, &path, compress)?;
println!(" Sen's slope → {}", path.display());
}
_ => anyhow::bail!(
"unknown trend method: '{}'. Use: linear, mann-kendall, sens",
method
),
}
done("Trend analysis", &outdir, start.elapsed());
Ok(())
}
fn handle_change(
before: PathBuf,
after: PathBuf,
output: PathBuf,
decrease_threshold: f64,
increase_threshold: f64,
compress: bool,
) -> Result<()> {
let before_r = read_dem(&before)?;
let after_r = read_dem(&after)?;
let start = Instant::now();
println!("Computing change detection...");
let params = RasterDiffParams {
decrease_threshold,
increase_threshold,
};
let (diff, cat) =
raster_difference(&before_r, &after_r, params).context("raster_difference failed")?;
write_result(&diff, &output, compress)?;
let cat_path = output.with_file_name(format!(
"{}_categorical.tif",
output.file_stem().unwrap_or_default().to_string_lossy()
));
write_result(&cat, &cat_path, compress)?;
println!(" difference → {}", output.display());
println!(
" categorical → {} (1=decrease, 2=stable, 3=increase)",
cat_path.display()
);
done("Change detection", &output, start.elapsed());
Ok(())
}
fn handle_anomaly(
reference: Vec<PathBuf>,
target: Vec<PathBuf>,
outdir: PathBuf,
method_str: String,
compress: bool,
) -> Result<()> {
let ref_rasters = load_raster_stack(&reference)?;
let tgt_rasters = load_raster_stack(&target)?;
let ref_refs: Vec<&surtgis_core::Raster<f64>> = ref_rasters.iter().collect();
let tgt_refs: Vec<&surtgis_core::Raster<f64>> = tgt_rasters.iter().collect();
let method = AnomalyMethod::from_str(&method_str).map_err(|e| anyhow::anyhow!("{}", e))?;
std::fs::create_dir_all(&outdir)?;
let start = Instant::now();
println!(
"Computing anomaly ({:?}) with {} reference, {} target rasters...",
method,
ref_rasters.len(),
tgt_rasters.len()
);
let results =
temporal_anomaly(&ref_refs, &tgt_refs, method).map_err(|e| anyhow::anyhow!("{}", e))?;
for (i, result) in results.iter().enumerate() {
let name = target[i]
.file_stem()
.map(|s| s.to_string_lossy().to_string())
.unwrap_or_else(|| format!("anomaly_{}", i));
let path = outdir.join(format!("{}_anomaly.tif", name));
write_result(result, &path, compress)?;
println!(" {} → {}", name, path.display());
}
done("Anomaly detection", &outdir, start.elapsed());
Ok(())
}
fn handle_phenology(
input: Vec<PathBuf>,
outdir: PathBuf,
doys_str: Option<String>,
threshold: f64,
smooth: usize,
compress: bool,
) -> Result<()> {
let rasters = load_raster_stack(&input)?;
let refs: Vec<&surtgis_core::Raster<f64>> = rasters.iter().collect();
let doys = doys_str.as_ref().map(|s| parse_times(s)).transpose()?;
let params = PhenologyParams {
threshold,
smooth_window: if smooth % 2 == 1 { smooth } else { smooth + 1 },
};
std::fs::create_dir_all(&outdir)?;
let start = Instant::now();
println!(
"Extracting vegetation phenology ({} time steps, threshold={:.2})...",
rasters.len(),
threshold
);
let result = vegetation_phenology(&refs, doys.as_deref(), ¶ms)
.map_err(|e| anyhow::anyhow!("{}", e))?;
let outputs = [
(&result.sos, "sos", "Start of Season"),
(&result.eos, "eos", "End of Season"),
(&result.peak, "peak", "Peak value"),
(&result.peak_time, "peak_time", "Peak time"),
(&result.amplitude, "amplitude", "Amplitude"),
(&result.season_length, "season_length", "Season length"),
];
for (raster, name, desc) in &outputs {
let path = outdir.join(format!("{}.tif", name));
write_result(raster, &path, compress)?;
println!(" {} ({}) → {}", name, desc, path.display());
}
done("Phenology", &outdir, start.elapsed());
Ok(())
}