use std::path::Path;
use std::time::Instant;
use anyhow::{Context, Result};
use geo::Point;
use surtgis_algorithms::interpolation::{
DriftOrder, FittedVariogram, OrdinaryKrigingParams, RegressionKrigingParams, SamplePoint,
UniversalKrigingParams, VariogramModel, VariogramParams, empirical_variogram,
fit_best_variogram, ordinary_kriging, regression_kriging, universal_kriging,
};
use surtgis_core::io::read_geotiff;
use surtgis_core::raster::Raster;
use surtgis_core::vector::{AttributeValue, read_vector};
use crate::commands::InterpolationCommands;
use crate::helpers::{done, spinner, write_result};
pub fn handle(action: InterpolationCommands, compress: bool) -> Result<()> {
match action {
InterpolationCommands::Variogram {
points,
attribute,
bins,
max_lag,
output,
} => {
let start = Instant::now();
let samples = load_samples(&points, &attribute)?;
println!("Loaded {} sample points", samples.len());
let pb = spinner("Computing variogram...");
let params = VariogramParams {
n_lags: bins,
max_lag,
lag_tolerance: 1.0,
};
let empirical = empirical_variogram(&samples, params)
.context("Failed to compute empirical variogram")?;
let fitted = fit_best_variogram(&empirical).context("Failed to fit variogram model")?;
pb.finish_and_clear();
println!(
"Model: {:?} — range={:.2}, sill={:.6}, nugget={:.6}",
fitted.model, fitted.range, fitted.sill, fitted.nugget
);
let result = serde_json::json!({
"empirical": {
"lags": empirical.lags,
"semivariance": empirical.semivariance,
"pair_counts": empirical.pair_counts,
},
"fitted": {
"model": format!("{:?}", fitted.model),
"range": fitted.range,
"sill": fitted.sill,
"nugget": fitted.nugget,
"partial_sill": fitted.partial_sill,
},
"n_points": samples.len(),
});
std::fs::write(&output, serde_json::to_string_pretty(&result)?)
.context("Failed to write variogram JSON")?;
done("Variogram", &output, start.elapsed());
}
InterpolationCommands::Kriging {
points,
attribute,
reference,
model,
range,
sill,
nugget,
max_neighbors,
output,
} => {
let start = Instant::now();
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let variogram = if let (Some(r), Some(s)) = (range, sill) {
FittedVariogram {
model: parse_variogram_model(&model),
range: r,
sill: s,
nugget,
partial_sill: s - nugget,
rss: 0.0,
}
} else {
let pb = spinner("Auto-fitting variogram...");
let emp = empirical_variogram(
&samples,
VariogramParams {
n_lags: 15,
max_lag: None,
lag_tolerance: 1.0,
},
)?;
let v = fit_best_variogram(&emp)?;
pb.finish_and_clear();
println!(
"Auto-fit: {:?} (range={:.1}, sill={:.4})",
v.model, v.range, v.sill
);
v
};
let pb = spinner("Kriging...");
let params = OrdinaryKrigingParams {
rows,
cols,
transform: ref_raster.transform().clone(),
max_points: max_neighbors,
max_radius: None,
compute_variance: false,
};
let result =
ordinary_kriging(&samples, &variogram, params).context("Kriging failed")?;
pb.finish_and_clear();
let mut out = result.estimate;
out.set_transform(ref_raster.transform().clone());
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("Ordinary Kriging", &output, start.elapsed());
}
InterpolationCommands::UniversalKriging {
points,
attribute,
reference,
drift,
model,
max_neighbors,
output,
} => {
let start = Instant::now();
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
let drift_order = match drift.to_lowercase().as_str() {
"quadratic" => DriftOrder::Quadratic,
_ => DriftOrder::Linear,
};
let pb = spinner("Auto-fitting variogram...");
let emp = empirical_variogram(
&samples,
VariogramParams {
n_lags: 15,
max_lag: None,
lag_tolerance: 1.0,
},
)?;
let variogram = fit_best_variogram(&emp)?;
pb.finish_and_clear();
let pb = spinner("Universal Kriging...");
let params = UniversalKrigingParams {
rows,
cols,
transform: ref_raster.transform().clone(),
max_points: max_neighbors,
max_radius: None,
drift_order,
compute_variance: false,
};
let result = universal_kriging(&samples, &variogram, params)
.context("Universal Kriging failed")?;
pb.finish_and_clear();
let mut out = result.estimate;
out.set_transform(ref_raster.transform().clone());
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("Universal Kriging", &output, start.elapsed());
}
InterpolationCommands::RegressionKriging {
points,
attribute,
features: _features,
reference,
model: _model,
output,
} => {
let start = Instant::now();
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let pb = spinner("Regression-Kriging...");
let params = RegressionKrigingParams {
rows,
cols,
transform: ref_raster.transform().clone(),
max_points: 16,
max_radius: None,
compute_variance: false,
variogram_params: VariogramParams {
n_lags: 15,
max_lag: None,
lag_tolerance: 1.0,
},
};
let result =
regression_kriging(&samples, params).context("Regression-Kriging failed")?;
pb.finish_and_clear();
let mut out = result.estimate;
out.set_transform(ref_raster.transform().clone());
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("Regression-Kriging", &output, start.elapsed());
}
InterpolationCommands::Idw {
points,
attribute,
reference,
power,
max_radius,
max_points,
output,
} => {
use surtgis_algorithms::interpolation::{IdwParams, idw};
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let start = Instant::now();
let pb = spinner(&format!("IDW (power={})...", power));
let result = idw(
&samples,
IdwParams {
power,
max_radius,
max_points,
rows,
cols,
transform: ref_raster.transform().clone(),
..IdwParams::default()
},
)
.context("IDW failed")?;
pb.finish_and_clear();
let mut out = result;
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("IDW", &output, start.elapsed());
}
InterpolationCommands::NearestNeighbor {
points,
attribute,
reference,
max_radius,
output,
} => {
use surtgis_algorithms::interpolation::{NearestNeighborParams, nearest_neighbor};
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let start = Instant::now();
let pb = spinner("Nearest Neighbor...");
let result = nearest_neighbor(
&samples,
NearestNeighborParams {
max_radius,
rows,
cols,
transform: ref_raster.transform().clone(),
},
)
.context("Nearest Neighbor failed")?;
pb.finish_and_clear();
let mut out = result;
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("Nearest Neighbor", &output, start.elapsed());
}
InterpolationCommands::NaturalNeighbor {
points,
attribute,
reference,
output,
} => {
use surtgis_algorithms::interpolation::{NaturalNeighborParams, natural_neighbor};
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let start = Instant::now();
let pb = spinner("Natural Neighbor...");
let result = natural_neighbor(
&samples,
NaturalNeighborParams {
max_neighbors: 20,
sub_resolution: 11,
rows,
cols,
transform: ref_raster.transform().clone(),
},
)
.context("Natural Neighbor failed")?;
pb.finish_and_clear();
let mut out = result;
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("Natural Neighbor", &output, start.elapsed());
}
InterpolationCommands::Tps {
points,
attribute,
reference,
smoothing,
output,
} => {
use surtgis_algorithms::interpolation::{TpsParams, tps_interpolation};
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let start = Instant::now();
let pb = spinner("Thin Plate Spline...");
let result = tps_interpolation(
&samples,
TpsParams {
rows,
cols,
transform: ref_raster.transform().clone(),
smoothing,
},
)
.context("TPS failed")?;
pb.finish_and_clear();
let mut out = result;
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("TPS", &output, start.elapsed());
}
InterpolationCommands::Tin {
points,
attribute,
reference,
output,
} => {
use surtgis_algorithms::interpolation::{TinParams, tin_interpolation};
let samples = load_samples(&points, &attribute)?;
let ref_raster: Raster<f64> = read_geotiff(&reference, None)?;
let (rows, cols) = ref_raster.shape();
println!("Points: {}, Grid: {}x{}", samples.len(), cols, rows);
let start = Instant::now();
let pb = spinner("TIN interpolation...");
let result = tin_interpolation(
&samples,
TinParams {
rows,
cols,
transform: ref_raster.transform().clone(),
},
)
.context("TIN failed")?;
pb.finish_and_clear();
let mut out = result;
out.set_crs(ref_raster.crs().cloned());
write_result(&out, &output, compress)?;
done("TIN", &output, start.elapsed());
}
}
Ok(())
}
fn load_samples(path: &Path, attribute: &str) -> Result<Vec<SamplePoint>> {
let fc = read_vector(path).context("Failed to read vector file")?;
let mut samples = Vec::new();
for feature in &fc.features {
let geom = feature
.geometry
.as_ref()
.ok_or_else(|| anyhow::anyhow!("Feature missing geometry"))?;
let (x, y) = match geom {
geo::Geometry::Point(p) => (p.x(), p.y()),
_ => continue, };
let value = match feature.properties.get(attribute) {
Some(AttributeValue::Float(v)) => *v,
Some(AttributeValue::Int(v)) => *v as f64,
_ => anyhow::bail!("Missing or non-numeric attribute '{}'", attribute),
};
samples.push(SamplePoint { x, y, value });
}
anyhow::ensure!(!samples.is_empty(), "No valid point samples found");
Ok(samples)
}
fn parse_variogram_model(s: &str) -> VariogramModel {
match s.to_lowercase().as_str() {
"exponential" => VariogramModel::Exponential,
"gaussian" => VariogramModel::Gaussian,
_ => VariogramModel::Spherical,
}
}