use anyhow::{Context, Result};
use std::time::Instant;
use surtgis_algorithms::terrain::{
advanced_curvatures, aspect, convergence_index, curvature, dev, eastness, geomorphons,
hillshade, landform_classification, ls_factor, mrvbf, multidirectional_hillshade,
negative_openness, northness, positive_openness, relative_slope_position, sky_view_factor,
slope, surface_area_ratio, tpi, tri, valley_depth, viewshed, vrm,
AspectOutput, AspectStreaming, ConvergenceParams, ConvergenceStreaming,
CurvatureParams, CurvatureStreaming, CurvatureType, DevParams, DevStreaming,
EastnessStreaming, GeomorphonParams, HillshadeParams, HillshadeStreaming, LandformParams,
LsFactorParams, MultiHillshadeParams, MultiHillshadeStreaming, MrvbfParams,
NorthnessStreaming, OpennessParams, SarParams, SlopeParams, SlopeStreaming, SlopeUnits,
SvfParams, TpiParams, TpiStreaming, TriParams, TriStreaming, ViewshedParams, VrmParams,
VrmStreaming,
};
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::{solar_radiation, SolarParams};
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::{solar_radiation_annual, SolarParams};
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(&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.annual.total, &output, compress)?;
done("Annual solar radiation", &output, start.elapsed());
}
TerrainCommands::ContourLines { input, output, interval, base } => {
use surtgis_algorithms::terrain::{contour_lines, ContourParams};
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::{cost_distance, CostDistanceParams};
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::{gaussian_smoothing, GaussianSmoothingParams};
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::{feature_preserving_smoothing, SmoothingParams};
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::{wind_exposure, WindExposureParams};
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::{sti, StiParams};
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::{uncertainty, UncertaintyParams};
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::{viewshed_pderl, PderlViewshedParams};
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::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(())
}