use anyhow::{Context, Result};
use std::time::Instant;
use surtgis_algorithms::terrain::{
AspectOutput, AspectStreaming, CircularVarianceParams, CircularVarianceStreaming,
ConvergenceParams, ConvergenceStreaming, CurvatureParams, CurvatureStreaming, CurvatureType,
DevParams, DevStreaming, DiffFromMeanParams, DiffFromMeanStreaming, DirectionalReliefParams,
DownslopeIndexParams, EastnessStreaming, EdgeDensityParams, GeomorphonParams, HillshadeParams,
HillshadeStreaming, LandformParams, LsFactorParams, MrvbfParams, MultiHillshadeParams,
MultiHillshadeStreaming, NormalDeviationParams, NormalDeviationStreaming, NorthnessStreaming,
OpennessParams, PennockParams, PercentElevRangeParams, PercentElevRangeStreaming,
RelativeAspectParams, SarParams, SlopeParams, SlopeStreaming, SlopeUnits,
SphericalStdDevParams, SphericalStdDevStreaming, SvfParams, TpiParams, TpiStreaming, TriParams,
TriStreaming, ViewshedParams, VrmParams, VrmStreaming, advanced_curvatures, aspect,
circular_variance_aspect, convergence_index, curvature, dev, diff_from_mean_elev,
directional_relief, downslope_index, eastness, edge_density, elev_above_pit,
elev_relative_to_min_max, geomorphons, hillshade, hypsometric_hillshade,
landform_classification, ls_factor, max_branch_length, mrvbf, multidirectional_hillshade,
negative_openness, neighbour_stats, normal_vector_deviation, northness, pennock,
percent_elev_range, positive_openness, relative_aspect, relative_slope_position,
sky_view_factor, slope, spherical_std_dev, surface_area_ratio, tpi, tri, valley_depth,
viewshed, vrm,
};
use surtgis_core::StripProcessor;
use crate::commands::TerrainCommands;
use crate::helpers::{
done, parse_advanced_curvature_type, read_dem, spinner, write_result, write_result_u8,
};
use crate::memory;
pub fn handle(
algorithm: TerrainCommands,
compress: bool,
streaming: bool,
mem_limit_bytes: Option<u64>,
) -> Result<()> {
match algorithm {
TerrainCommands::Slope {
input,
output,
units,
z_factor,
} => {
let units = match units.to_lowercase().as_str() {
"degrees" | "deg" | "d" => SlopeUnits::Degrees,
"percent" | "pct" | "%" => SlopeUnits::Percent,
"radians" | "rad" | "r" => SlopeUnits::Radians,
_ => {
eprintln!("Unknown units: {}. Using degrees.", units);
SlopeUnits::Degrees
}
};
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = SlopeStreaming { units, z_factor };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to calculate slope (streaming)")?;
let elapsed = start.elapsed();
println!("Slope (streaming): {} x {}", cols, rows);
done("Slope", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = slope(&dem, SlopeParams { units, z_factor })
.context("Failed to calculate slope")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Slope", &output, elapsed);
}
}
TerrainCommands::Aspect {
input,
output,
format,
} => {
let fmt = match format.to_lowercase().as_str() {
"degrees" | "deg" | "d" => AspectOutput::Degrees,
"radians" | "rad" | "r" => AspectOutput::Radians,
"compass" | "c" => AspectOutput::Compass,
_ => {
eprintln!("Unknown format: {}. Using degrees.", format);
AspectOutput::Degrees
}
};
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = AspectStreaming { output_format: fmt };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to calculate aspect (streaming)")?;
let elapsed = start.elapsed();
println!("Aspect (streaming): {} x {}", cols, rows);
done("Aspect", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = aspect(&dem, fmt).context("Failed to calculate aspect")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Aspect", &output, elapsed);
}
}
TerrainCommands::Hillshade {
input,
output,
azimuth,
altitude,
z_factor,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = HillshadeStreaming {
azimuth,
altitude,
z_factor,
};
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to calculate hillshade (streaming)")?;
let elapsed = start.elapsed();
println!("Hillshade (streaming): {} x {}", cols, rows);
done("Hillshade", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = hillshade(
&dem,
HillshadeParams {
azimuth,
altitude,
z_factor,
normalized: false,
},
)
.context("Failed to calculate hillshade")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Hillshade", &output, elapsed);
}
}
TerrainCommands::Curvature {
input,
output,
curvature_type,
z_factor,
} => {
let ct = match curvature_type.to_lowercase().as_str() {
"general" | "mean" | "g" => CurvatureType::General,
"profile" | "prof" | "p" => CurvatureType::Profile,
"plan" | "tangential" | "t" => CurvatureType::Plan,
_ => {
eprintln!("Unknown type: {}. Using general.", curvature_type);
CurvatureType::General
}
};
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = CurvatureStreaming {
curvature_type: ct,
z_factor,
..Default::default()
};
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to calculate curvature (streaming)")?;
let elapsed = start.elapsed();
println!("Curvature (streaming): {} x {}", cols, rows);
done("Curvature", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = curvature(
&dem,
CurvatureParams {
curvature_type: ct,
z_factor,
..Default::default()
},
)
.context("Failed to calculate curvature")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Curvature", &output, elapsed);
}
}
TerrainCommands::Tpi {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = TpiStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to calculate TPI (streaming)")?;
let elapsed = start.elapsed();
println!("TPI (streaming): {} x {}", cols, rows);
done("TPI", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = tpi(&dem, TpiParams { radius }).context("Failed to calculate TPI")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("TPI", &output, elapsed);
}
}
TerrainCommands::Tri {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = TriStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to calculate TRI (streaming)")?;
let elapsed = start.elapsed();
println!("TRI (streaming): {} x {}", cols, rows);
done("TRI", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = tri(&dem, TriParams { radius }).context("Failed to calculate TRI")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("TRI", &output, elapsed);
}
}
TerrainCommands::Landform {
input,
output,
small_radius,
large_radius,
threshold,
slope_threshold,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = landform_classification(
&dem,
LandformParams {
small_radius,
large_radius,
tpi_threshold: threshold,
slope_threshold,
},
)
.context("Failed to classify landforms")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Landform classification", &output, elapsed);
}
TerrainCommands::Geomorphons {
input,
output,
radius,
flatness,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = geomorphons(
&dem,
GeomorphonParams {
radius,
flatness_threshold: flatness,
},
)
.context("Failed to compute geomorphons")?;
let elapsed = start.elapsed();
write_result_u8(&result, &output, compress)?;
done("Geomorphons", &output, elapsed);
}
TerrainCommands::Northness { input, output } => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = NorthnessStreaming;
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute northness (streaming)")?;
let elapsed = start.elapsed();
println!("Northness (streaming): {} x {}", cols, rows);
done("Northness", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = northness(&dem).context("Failed to compute northness")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Northness", &output, elapsed);
}
}
TerrainCommands::Eastness { input, output } => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = EastnessStreaming;
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute eastness (streaming)")?;
let elapsed = start.elapsed();
println!("Eastness (streaming): {} x {}", cols, rows);
done("Eastness", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = eastness(&dem).context("Failed to compute eastness")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Eastness", &output, elapsed);
}
}
TerrainCommands::OpennessPositive {
input,
output,
radius,
directions,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = positive_openness(&dem, OpennessParams { radius, directions })
.context("Failed to compute positive openness")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Positive openness", &output, elapsed);
}
TerrainCommands::OpennessNegative {
input,
output,
radius,
directions,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = negative_openness(&dem, OpennessParams { radius, directions })
.context("Failed to compute negative openness")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Negative openness", &output, elapsed);
}
TerrainCommands::Svf {
input,
output,
radius,
directions,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = sky_view_factor(&dem, SvfParams { radius, directions })
.context("Failed to compute SVF")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Sky View Factor", &output, elapsed);
}
TerrainCommands::Mrvbf {
input,
output,
mrrtf_output,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let (mrvbf_r, mrrtf_r) =
mrvbf(&dem, MrvbfParams::default()).context("Failed to compute MRVBF")?;
let elapsed = start.elapsed();
write_result(&mrvbf_r, &output, compress)?;
done("MRVBF", &output, elapsed);
if let Some(mrrtf_path) = mrrtf_output {
write_result(&mrrtf_r, &mrrtf_path, compress)?;
println!("MRRTF saved to: {}", mrrtf_path.display());
}
}
TerrainCommands::Dev {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = DevStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute DEV (streaming)")?;
let elapsed = start.elapsed();
println!("DEV (streaming): {} x {}", cols, rows);
done("DEV", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = dev(&dem, DevParams { radius }).context("Failed to compute DEV")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("DEV", &output, elapsed);
}
}
TerrainCommands::Vrm {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = VrmStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute VRM (streaming)")?;
let elapsed = start.elapsed();
println!("VRM (streaming): {} x {}", cols, rows);
done("VRM", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = vrm(&dem, VrmParams { radius }).context("Failed to compute VRM")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("VRM", &output, elapsed);
}
}
TerrainCommands::AdvancedCurvature {
input,
output,
curvature_type,
} => {
let ct = parse_advanced_curvature_type(&curvature_type)?;
let dem = read_dem(&input)?;
let start = Instant::now();
let result =
advanced_curvatures(&dem, ct).context("Failed to compute advanced curvature")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Advanced curvature", &output, elapsed);
}
TerrainCommands::Viewshed {
input,
output,
observer_row,
observer_col,
observer_height,
target_height,
max_radius,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = viewshed(
&dem,
ViewshedParams {
observer_row,
observer_col,
observer_height,
target_height,
max_radius,
},
)
.context("Failed to compute viewshed")?;
let elapsed = start.elapsed();
write_result_u8(&result, &output, compress)?;
done("Viewshed", &output, elapsed);
}
TerrainCommands::Convergence {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = ConvergenceStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute convergence index (streaming)")?;
let elapsed = start.elapsed();
println!("Convergence (streaming): {} x {}", cols, rows);
done("Convergence index", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = convergence_index(&dem, ConvergenceParams { radius })
.context("Failed to compute convergence index")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Convergence index", &output, elapsed);
}
}
TerrainCommands::MultiHillshade { input, output } => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = MultiHillshadeStreaming::default();
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute multi-directional hillshade (streaming)")?;
let elapsed = start.elapsed();
println!(
"Multi-directional hillshade (streaming): {} x {}",
cols, rows
);
done("Multi-directional hillshade", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = multidirectional_hillshade(&dem, MultiHillshadeParams::default())
.context("Failed to compute multi-directional hillshade")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Multi-directional hillshade", &output, elapsed);
}
}
TerrainCommands::LsFactor {
flow_acc,
slope: slope_path,
output,
cell_size,
} => {
let facc = read_dem(&flow_acc)?;
let slp = read_dem(&slope_path)?;
let start = Instant::now();
let result = ls_factor(
&facc,
&slp,
LsFactorParams {
cell_size,
..Default::default()
},
)
.context("Failed to compute LS-Factor")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("LS-Factor", &output, elapsed);
}
TerrainCommands::ValleyDepth { input, output } => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = valley_depth(&dem).context("Failed to compute valley depth")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Valley depth", &output, elapsed);
}
TerrainCommands::RelativeSlopePosition {
hand,
valley_depth: vd_path,
output,
} => {
let hand_r = read_dem(&hand)?;
let vd_r = read_dem(&vd_path)?;
let start = Instant::now();
let result = relative_slope_position(&hand_r, &vd_r)
.context("Failed to compute relative slope position")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Relative slope position", &output, elapsed);
}
TerrainCommands::SurfaceAreaRatio {
input,
output,
radius,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = surface_area_ratio(&dem, SarParams { radius })
.context("Failed to compute surface area ratio")?;
let elapsed = start.elapsed();
write_result(&result, &output, compress)?;
done("Surface area ratio", &output, elapsed);
}
TerrainCommands::SolarRadiation {
input,
output,
day,
hour: _hour,
latitude,
} => {
use surtgis_algorithms::terrain::{SolarParams, solar_radiation};
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("Computing slope+aspect for solar...");
let slope_r = slope(
&dem,
SlopeParams {
units: SlopeUnits::Radians,
z_factor: 1.0,
},
)
.context("Slope failed")?;
let aspect_r = aspect(&dem, AspectOutput::Radians).context("Aspect failed")?;
pb.finish_and_clear();
let pb = spinner(&format!("Solar radiation (day={})...", day));
let result = solar_radiation(
&slope_r,
&aspect_r,
SolarParams {
day,
latitude,
transmittance: 0.7,
..SolarParams::default()
},
)
.context("Solar radiation failed")?;
pb.finish_and_clear();
write_result(&result.total, &output, compress)?;
done("Solar radiation", &output, start.elapsed());
}
TerrainCommands::SolarRadiationAnnual {
input,
output,
latitude,
} => {
use surtgis_algorithms::terrain::{SolarParams, solar_radiation_annual_only};
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("Computing slope+aspect...");
let slope_r = slope(
&dem,
SlopeParams {
units: SlopeUnits::Radians,
z_factor: 1.0,
},
)
.context("Slope failed")?;
let aspect_r = aspect(&dem, AspectOutput::Radians).context("Aspect failed")?;
pb.finish_and_clear();
let pb = spinner("Annual solar radiation...");
let result = solar_radiation_annual_only(
&slope_r,
&aspect_r,
SolarParams {
day: 1,
latitude,
transmittance: 0.7,
..SolarParams::default()
},
)
.context("Annual solar radiation failed")?;
pb.finish_and_clear();
write_result(&result.total, &output, compress)?;
done("Annual solar radiation", &output, start.elapsed());
}
TerrainCommands::ContourLines {
input,
output,
interval,
base,
} => {
use surtgis_algorithms::terrain::{ContourParams, contour_lines};
let dem = read_dem(&input)?;
let start = Instant::now();
let result = contour_lines(&dem, ContourParams { interval, base })
.context("Contour lines failed")?;
write_result(&result, &output, compress)?;
done("Contour lines", &output, start.elapsed());
}
TerrainCommands::CostDistance {
input,
sources,
output,
} => {
use surtgis_algorithms::terrain::{CostDistanceParams, cost_distance};
let cost = read_dem(&input)?;
let pb = spinner("Reading sources...");
let src: surtgis_core::Raster<f64> = surtgis_core::io::read_geotiff(&sources, None)
.context("Failed to read source raster")?;
pb.finish_and_clear();
let mut src_positions = Vec::new();
let (sr, sc) = src.shape();
for r in 0..sr {
for c in 0..sc {
if let Some(&v) = src.data().get([r, c]) {
if v != 0.0 && v.is_finite() {
src_positions.push((r, c));
}
}
}
}
println!(" {} source cells", src_positions.len());
let start = Instant::now();
let pb = spinner("Cost distance...");
let result = cost_distance(
&cost,
CostDistanceParams {
sources: src_positions,
},
)
.context("Cost distance failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Cost distance", &output, start.elapsed());
}
TerrainCommands::ShapeIndex { input, output } => {
use surtgis_algorithms::terrain::shape_index;
let dem = read_dem(&input)?;
let start = Instant::now();
let result = shape_index(&dem).context("Shape index failed")?;
write_result(&result, &output, compress)?;
done("Shape index", &output, start.elapsed());
}
TerrainCommands::Curvedness { input, output } => {
use surtgis_algorithms::terrain::curvedness;
let dem = read_dem(&input)?;
let start = Instant::now();
let result = curvedness(&dem).context("Curvedness failed")?;
write_result(&result, &output, compress)?;
done("Curvedness", &output, start.elapsed());
}
TerrainCommands::GaussianSmoothing {
input,
output,
sigma,
radius,
} => {
use surtgis_algorithms::terrain::{GaussianSmoothingParams, gaussian_smoothing};
let dem = read_dem(&input)?;
let r = radius.unwrap_or((3.0 * sigma).ceil() as usize);
let start = Instant::now();
let pb = spinner(&format!(
"Gaussian smoothing (sigma={}, radius={})...",
sigma, r
));
let result = gaussian_smoothing(&dem, GaussianSmoothingParams { radius: r, sigma })
.context("Smoothing failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Gaussian smoothing", &output, start.elapsed());
}
TerrainCommands::FeaturePreservingSmoothing {
input,
output,
strength: _strength,
iterations,
} => {
use surtgis_algorithms::terrain::{SmoothingParams, feature_preserving_smoothing};
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("Feature-preserving smoothing...");
let result = feature_preserving_smoothing(
&dem,
SmoothingParams {
radius: 2,
iterations,
threshold: 15.0,
},
)
.context("Smoothing failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Feature-preserving smoothing", &output, start.elapsed());
}
TerrainCommands::WindExposure {
input,
output,
direction,
radius,
} => {
use surtgis_algorithms::terrain::{WindExposureParams, wind_exposure};
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner(&format!(
"Wind exposure (dir={}°, radius={})...",
direction, radius
));
let result = wind_exposure(
&dem,
WindExposureParams {
radius,
directions: 8,
wind_direction: Some(direction.to_radians()),
wind_window: 45.0,
},
)
.context("Wind exposure failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Wind exposure", &output, start.elapsed());
}
TerrainCommands::HorizonAngle {
input,
output,
azimuth,
radius,
} => {
use surtgis_algorithms::terrain::horizon_angle_map;
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner(&format!(
"Horizon angle (az={}°, radius={})...",
azimuth, radius
));
let result = horizon_angle_map(&dem, azimuth.to_radians(), radius)
.context("Horizon angle failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Horizon angle", &output, start.elapsed());
}
TerrainCommands::AccumulationZones { input, output } => {
use surtgis_algorithms::terrain::accumulation_zones;
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("Accumulation zones...");
let result = accumulation_zones(&dem).context("Accumulation zones failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Accumulation zones", &output, start.elapsed());
}
TerrainCommands::Spi {
flow_acc,
slope: slope_path,
output,
} => {
use surtgis_algorithms::terrain::spi;
let fa = read_dem(&flow_acc)?;
let sl = read_dem(&slope_path)?;
let start = Instant::now();
let result = spi(&fa, &sl).context("SPI failed")?;
write_result(&result, &output, compress)?;
done("SPI", &output, start.elapsed());
}
TerrainCommands::Sti {
flow_acc,
slope: slope_path,
output,
} => {
use surtgis_algorithms::terrain::{StiParams, sti};
let fa = read_dem(&flow_acc)?;
let sl = read_dem(&slope_path)?;
let start = Instant::now();
let result = sti(&fa, &sl, StiParams { m: 0.4, n: 1.3 }).context("STI failed")?;
write_result(&result, &output, compress)?;
done("STI", &output, start.elapsed());
}
TerrainCommands::Twi {
flow_acc,
slope: slope_path,
output,
} => {
use surtgis_algorithms::terrain::twi;
let fa = read_dem(&flow_acc)?;
let sl = read_dem(&slope_path)?;
let start = Instant::now();
let result = twi(&fa, &sl).context("TWI failed")?;
write_result(&result, &output, compress)?;
done("TWI", &output, start.elapsed());
}
TerrainCommands::LogTransform { input, output } => {
use surtgis_algorithms::terrain::log_transform;
let raster = read_dem(&input)?;
let start = Instant::now();
let result = log_transform(&raster).context("Log transform failed")?;
write_result(&result, &output, compress)?;
done("Log transform", &output, start.elapsed());
}
TerrainCommands::Uncertainty {
input,
outdir,
error_std,
n_simulations: _n,
} => {
use surtgis_algorithms::terrain::{UncertaintyParams, uncertainty};
std::fs::create_dir_all(&outdir).context("Failed to create output directory")?;
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("Uncertainty analysis...");
let result = uncertainty(
&dem,
UncertaintyParams {
dem_rmse: error_std,
},
)
.context("Uncertainty failed")?;
pb.finish_and_clear();
write_result(&result.slope_rmse, &outdir.join("slope_rmse.tif"), compress)?;
write_result(
&result.aspect_rmse,
&outdir.join("aspect_rmse.tif"),
compress,
)?;
done("Uncertainty", &outdir, start.elapsed());
}
TerrainCommands::ViewshedPderl {
input,
output,
row,
col,
height,
} => {
use surtgis_algorithms::terrain::{PderlViewshedParams, viewshed_pderl};
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("PDERL viewshed...");
let result = viewshed_pderl(
&dem,
PderlViewshedParams {
observer_row: row,
observer_col: col,
observer_height: height,
target_height: 0.0,
max_radius: 0,
},
)
.context("Viewshed failed")?;
pb.finish_and_clear();
write_result_u8(&result, &output, compress)?;
done("PDERL Viewshed", &output, start.elapsed());
}
TerrainCommands::ViewshedXdraw {
input,
output,
row,
col,
height,
} => {
use surtgis_algorithms::terrain::viewshed_xdraw;
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("XDraw viewshed...");
let result = viewshed_xdraw(
&dem,
ViewshedParams {
observer_row: row,
observer_col: col,
observer_height: height,
target_height: 0.0,
max_radius: 0,
},
)
.context("Viewshed failed")?;
pb.finish_and_clear();
write_result_u8(&result, &output, compress)?;
done("XDraw Viewshed", &output, start.elapsed());
}
TerrainCommands::ViewshedMultiple {
input,
output,
observers,
height,
} => {
use surtgis_algorithms::terrain::viewshed_multiple;
let dem = read_dem(&input)?;
let obs: Vec<(usize, usize)> = observers
.split(';')
.filter_map(|s| {
let parts: Vec<&str> = s.trim().split(',').collect();
if parts.len() == 2 {
Some((parts[0].trim().parse().ok()?, parts[1].trim().parse().ok()?))
} else {
None
}
})
.collect();
println!(" {} observer locations", obs.len());
let start = Instant::now();
let pb = spinner("Multiple viewshed...");
let result =
viewshed_multiple(&dem, &obs, height).context("Multiple viewshed failed")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Multiple Viewshed", &output, start.elapsed());
}
TerrainCommands::HypsometricHillshade {
input,
output,
azimuth,
altitude,
z_factor,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = hypsometric_hillshade(
&dem,
HillshadeParams {
azimuth,
altitude,
z_factor,
normalized: true,
},
)
.context("Failed to compute hypsometric hillshade")?;
write_result(&result, &output, compress)?;
done("Hypsometric Hillshade", &output, start.elapsed());
}
TerrainCommands::ElevRelative { input, output } => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = elev_relative_to_min_max(&dem)
.context("Failed to compute elevation relative to min/max")?;
write_result(&result, &output, compress)?;
done("Elev Relative to Min/Max", &output, start.elapsed());
}
TerrainCommands::DiffFromMean {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = DiffFromMeanStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute diff-from-mean (streaming)")?;
let elapsed = start.elapsed();
println!("Diff from Mean (streaming): {} x {}", cols, rows);
done("Diff from Mean Elev", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = diff_from_mean_elev(&dem, DiffFromMeanParams { radius })
.context("Failed to compute diff from mean elevation")?;
write_result(&result, &output, compress)?;
done("Diff from Mean Elev", &output, start.elapsed());
}
}
TerrainCommands::PercentElevRange {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = PercentElevRangeStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute percent-elev-range (streaming)")?;
let elapsed = start.elapsed();
println!("Percent Elev Range (streaming): {} x {}", cols, rows);
done("Percent Elev Range", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = percent_elev_range(&dem, PercentElevRangeParams { radius })
.context("Failed to compute percent elevation range")?;
write_result(&result, &output, compress)?;
done("Percent Elev Range", &output, start.elapsed());
}
}
TerrainCommands::ElevAbovePit { input, output } => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = elev_above_pit(&dem).context("Failed to compute elevation above pit")?;
write_result(&result, &output, compress)?;
done("Elev Above Pit", &output, start.elapsed());
}
TerrainCommands::CircularVarianceAspect {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = CircularVarianceStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute circular variance of aspect (streaming)")?;
let elapsed = start.elapsed();
println!("Circular Variance Aspect (streaming): {} x {}", cols, rows);
done("Circular Variance Aspect", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = circular_variance_aspect(&dem, CircularVarianceParams { radius })
.context("Failed to compute circular variance of aspect")?;
write_result(&result, &output, compress)?;
done("Circular Variance Aspect", &output, start.elapsed());
}
}
TerrainCommands::Neighbours {
input,
output,
stat,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = neighbour_stats(&dem).context("Failed to compute neighbour statistics")?;
if let Some(stat_name) = stat {
let raster = match stat_name.as_str() {
"max_downslope_change" => &result.max_downslope_change,
"min_downslope_change" => &result.min_downslope_change,
"max_upslope_change" => &result.max_upslope_change,
"num_downslope" => &result.num_downslope,
"num_upslope" => &result.num_upslope,
other => anyhow::bail!(
"Unknown stat '{}'. Valid: max_downslope_change, min_downslope_change, max_upslope_change, num_downslope, num_upslope",
other
),
};
write_result(raster, &output, compress)?;
done(
&format!("Neighbours ({})", stat_name),
&output,
start.elapsed(),
);
} else {
std::fs::create_dir_all(&output).context("Failed to create output directory")?;
write_result(
&result.max_downslope_change,
&output.join("max_downslope_change.tif"),
compress,
)?;
write_result(
&result.min_downslope_change,
&output.join("min_downslope_change.tif"),
compress,
)?;
write_result(
&result.max_upslope_change,
&output.join("max_upslope_change.tif"),
compress,
)?;
write_result(
&result.num_downslope,
&output.join("num_downslope.tif"),
compress,
)?;
write_result(
&result.num_upslope,
&output.join("num_upslope.tif"),
compress,
)?;
println!(" 5 neighbour statistics saved to {}", output.display());
done("Neighbour Stats", &output, start.elapsed());
}
}
TerrainCommands::Pennock {
input,
output,
slope_threshold,
curv_threshold,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = pennock(
&dem,
PennockParams {
slope_threshold,
profile_curv_threshold: curv_threshold,
plan_curv_threshold: 0.0,
},
)
.context("Failed to compute Pennock classification")?;
write_result(&result, &output, compress)?;
done("Pennock Landform", &output, start.elapsed());
}
TerrainCommands::EdgeDensity {
input,
output,
radius,
threshold,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = edge_density(&dem, EdgeDensityParams { radius, threshold })
.context("Failed to compute edge density")?;
write_result(&result, &output, compress)?;
done("Edge Density", &output, start.elapsed());
}
TerrainCommands::RelativeAspect {
input,
output,
sigma,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = relative_aspect(&dem, RelativeAspectParams { sigma })
.context("Failed to compute relative aspect")?;
write_result(&result, &output, compress)?;
done("Relative Aspect", &output, start.elapsed());
}
TerrainCommands::NormalDeviation {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = NormalDeviationStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute normal deviation (streaming)")?;
let elapsed = start.elapsed();
println!("Normal Deviation (streaming): {} x {}", cols, rows);
done("Normal Deviation", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = normal_vector_deviation(&dem, NormalDeviationParams { radius })
.context("Failed to compute normal vector deviation")?;
write_result(&result, &output, compress)?;
done("Normal Deviation", &output, start.elapsed());
}
}
TerrainCommands::SphericalStdDev {
input,
output,
radius,
} => {
let use_streaming = memory::should_stream(&input, mem_limit_bytes, streaming)?;
if use_streaming {
let algo = SphericalStdDevStreaming { radius };
let processor = StripProcessor::new(256);
let start = Instant::now();
let (rows, cols) = processor
.process(&input, &output, &algo, compress)
.context("Failed to compute spherical std dev (streaming)")?;
let elapsed = start.elapsed();
println!("Spherical Std Dev (streaming): {} x {}", cols, rows);
done("Spherical Std Dev", &output, elapsed);
} else {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = spherical_std_dev(&dem, SphericalStdDevParams { radius })
.context("Failed to compute spherical std dev of normals")?;
write_result(&result, &output, compress)?;
done("Spherical Std Dev", &output, start.elapsed());
}
}
TerrainCommands::DirectionalRelief {
input,
output,
radius,
azimuth,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = directional_relief(&dem, DirectionalReliefParams { radius, azimuth })
.context("Failed to compute directional relief")?;
write_result(&result, &output, compress)?;
done("Directional Relief", &output, start.elapsed());
}
TerrainCommands::DownslopeIndex {
input,
output,
drop: drop_val,
} => {
let dem = read_dem(&input)?;
let start = Instant::now();
let result = downslope_index(&dem, DownslopeIndexParams { drop: drop_val })
.context("Failed to compute downslope index")?;
write_result(&result, &output, compress)?;
done("Downslope Index", &output, start.elapsed());
}
TerrainCommands::MaxBranchLength { input, output } => {
let dem = read_dem(&input)?;
let start = Instant::now();
let pb = spinner("Computing max branch length (fill + flow dir + topo sort)...");
let result = max_branch_length(&dem).context("Failed to compute max branch length")?;
pb.finish_and_clear();
write_result(&result, &output, compress)?;
done("Max Branch Length", &output, start.elapsed());
}
TerrainCommands::All { input, outdir } => {
std::fs::create_dir_all(&outdir).context("Failed to create output directory")?;
let dem = read_dem(&input)?;
let start = Instant::now();
println!("Computing all terrain factors...");
let s = slope(
&dem,
SlopeParams {
units: SlopeUnits::Degrees,
z_factor: 1.0,
},
)
.context("slope")?;
write_result(&s, &outdir.join("slope.tif"), compress)?;
println!(" slope.tif");
let a = aspect(&dem, AspectOutput::Degrees).context("aspect")?;
write_result(&a, &outdir.join("aspect.tif"), compress)?;
println!(" aspect.tif");
let h = hillshade(
&dem,
HillshadeParams {
azimuth: 315.0,
altitude: 45.0,
z_factor: 1.0,
normalized: false,
},
)
.context("hillshade")?;
write_result(&h, &outdir.join("hillshade.tif"), compress)?;
println!(" hillshade.tif");
let n = northness(&dem).context("northness")?;
write_result(&n, &outdir.join("northness.tif"), compress)?;
println!(" northness.tif");
let e = eastness(&dem).context("eastness")?;
write_result(&e, &outdir.join("eastness.tif"), compress)?;
println!(" eastness.tif");
let c = curvature(
&dem,
CurvatureParams {
curvature_type: CurvatureType::General,
z_factor: 1.0,
..Default::default()
},
)
.context("curvature")?;
write_result(&c, &outdir.join("curvature.tif"), compress)?;
println!(" curvature.tif");
let t = tpi(&dem, TpiParams { radius: 10 }).context("tpi")?;
write_result(&t, &outdir.join("tpi.tif"), compress)?;
println!(" tpi.tif");
let tr = tri(&dem, TriParams { radius: 1 }).context("tri")?;
write_result(&tr, &outdir.join("tri.tif"), compress)?;
println!(" tri.tif");
let g = geomorphons(
&dem,
GeomorphonParams {
radius: 10,
flatness_threshold: 1.0,
},
)
.context("geomorphons")?;
write_result_u8(&g, &outdir.join("geomorphons.tif"), compress)?;
println!(" geomorphons.tif");
let d = dev(&dem, DevParams { radius: 10 }).context("dev")?;
write_result(&d, &outdir.join("dev.tif"), compress)?;
println!(" dev.tif");
let v = vrm(&dem, VrmParams { radius: 1 }).context("vrm")?;
write_result(&v, &outdir.join("vrm.tif"), compress)?;
println!(" vrm.tif");
let ci =
convergence_index(&dem, ConvergenceParams { radius: 3 }).context("convergence")?;
write_result(&ci, &outdir.join("convergence.tif"), compress)?;
println!(" convergence.tif");
let op = positive_openness(
&dem,
OpennessParams {
radius: 10,
directions: 8,
},
)
.context("openness_positive")?;
write_result(&op, &outdir.join("openness_positive.tif"), compress)?;
println!(" openness_positive.tif");
let on = negative_openness(
&dem,
OpennessParams {
radius: 10,
directions: 8,
},
)
.context("openness_negative")?;
write_result(&on, &outdir.join("openness_negative.tif"), compress)?;
println!(" openness_negative.tif");
let svf_r = sky_view_factor(
&dem,
SvfParams {
radius: 10,
directions: 16,
},
)
.context("svf")?;
write_result(&svf_r, &outdir.join("svf.tif"), compress)?;
println!(" svf.tif");
let (mrvbf_r, mrrtf_r) = mrvbf(&dem, MrvbfParams::default()).context("mrvbf")?;
write_result(&mrvbf_r, &outdir.join("mrvbf.tif"), compress)?;
write_result(&mrrtf_r, &outdir.join("mrrtf.tif"), compress)?;
println!(" mrvbf.tif, mrrtf.tif");
let elapsed = start.elapsed();
println!("\nAll terrain factors saved to: {}", outdir.display());
println!(" 17 products, processing time: {:.2?}", elapsed);
}
}
Ok(())
}
pub fn handle_terrain_all(
input: &std::path::Path,
outdir: &std::path::Path,
compress: bool,
) -> Result<()> {
std::fs::create_dir_all(outdir).context("Failed to create output directory")?;
let dem = read_dem(&input.to_path_buf())?;
let start = std::time::Instant::now();
println!("Computing all terrain factors...");
let s = slope(
&dem,
SlopeParams {
units: SlopeUnits::Degrees,
z_factor: 1.0,
},
)
.context("slope")?;
write_result(&s, &outdir.join("slope.tif"), compress)?;
println!(" slope.tif");
let a = aspect(&dem, AspectOutput::Degrees).context("aspect")?;
write_result(&a, &outdir.join("aspect.tif"), compress)?;
println!(" aspect.tif");
let h = hillshade(
&dem,
HillshadeParams {
azimuth: 315.0,
altitude: 45.0,
z_factor: 1.0,
normalized: false,
},
)
.context("hillshade")?;
write_result(&h, &outdir.join("hillshade.tif"), compress)?;
println!(" hillshade.tif");
let n = northness(&dem).context("northness")?;
write_result(&n, &outdir.join("northness.tif"), compress)?;
println!(" northness.tif");
let e = eastness(&dem).context("eastness")?;
write_result(&e, &outdir.join("eastness.tif"), compress)?;
println!(" eastness.tif");
let c = curvature(
&dem,
CurvatureParams {
curvature_type: CurvatureType::General,
z_factor: 1.0,
..Default::default()
},
)
.context("curvature")?;
write_result(&c, &outdir.join("curvature.tif"), compress)?;
println!(" curvature.tif");
let t = tpi(&dem, TpiParams { radius: 10 }).context("tpi")?;
write_result(&t, &outdir.join("tpi.tif"), compress)?;
println!(" tpi.tif");
let tr = tri(&dem, TriParams { radius: 1 }).context("tri")?;
write_result(&tr, &outdir.join("tri.tif"), compress)?;
println!(" tri.tif");
let g = geomorphons(
&dem,
GeomorphonParams {
radius: 10,
flatness_threshold: 1.0,
},
)
.context("geomorphons")?;
write_result_u8(&g, &outdir.join("geomorphons.tif"), compress)?;
println!(" geomorphons.tif");
let d = dev(&dem, DevParams { radius: 10 }).context("dev")?;
write_result(&d, &outdir.join("dev.tif"), compress)?;
println!(" dev.tif");
let v = vrm(&dem, VrmParams { radius: 1 }).context("vrm")?;
write_result(&v, &outdir.join("vrm.tif"), compress)?;
println!(" vrm.tif");
let ci = convergence_index(&dem, ConvergenceParams { radius: 3 }).context("convergence")?;
write_result(&ci, &outdir.join("convergence.tif"), compress)?;
println!(" convergence.tif");
let op = positive_openness(
&dem,
OpennessParams {
radius: 10,
directions: 8,
},
)
.context("openness_positive")?;
write_result(&op, &outdir.join("openness_positive.tif"), compress)?;
println!(" openness_positive.tif");
let on = negative_openness(
&dem,
OpennessParams {
radius: 10,
directions: 8,
},
)
.context("openness_negative")?;
write_result(&on, &outdir.join("openness_negative.tif"), compress)?;
println!(" openness_negative.tif");
let svf_r = sky_view_factor(
&dem,
SvfParams {
radius: 10,
directions: 16,
},
)
.context("svf")?;
write_result(&svf_r, &outdir.join("svf.tif"), compress)?;
println!(" svf.tif");
let (mrvbf_r, mrrtf_r) = mrvbf(&dem, MrvbfParams::default()).context("mrvbf")?;
write_result(&mrvbf_r, &outdir.join("mrvbf.tif"), compress)?;
write_result(&mrrtf_r, &outdir.join("mrrtf.tif"), compress)?;
println!(" mrvbf.tif, mrrtf.tif");
let elapsed = start.elapsed();
println!("\nAll terrain factors saved to: {}", outdir.display());
println!(" 17 products, processing time: {:.2?}", elapsed);
Ok(())
}