surtgis 0.6.24

High-performance geospatial analysis CLI
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

//! Handlers for classification algorithms.

use anyhow::{Context, Result};
use std::time::Instant;

use surtgis_algorithms::classification::{
    isodata, kmeans_raster, maximum_likelihood, minimum_distance, pca,
    signatures_from_training, IsodataParams, KmeansParams, PcaParams,
};
use surtgis_core::io::read_geotiff;

use crate::commands::ClassificationCommands;
use crate::helpers::{done, read_dem, spinner, write_result};

pub fn handle(algorithm: ClassificationCommands, compress: bool) -> Result<()> {
    match algorithm {
        ClassificationCommands::Kmeans {
            input, output, classes, max_iter, convergence, seed,
        } => {
            let raster = read_dem(&input)?;
            let start = Instant::now();
            let pb = spinner("K-means clustering...");
            let result = kmeans_raster(
                &raster,
                KmeansParams {
                    k: classes,
                    max_iterations: max_iter,
                    convergence,
                    seed,
                },
            )
            .context("K-means failed")?;
            pb.finish_and_clear();
            write_result(&result, &output, compress)?;
            done("K-means", &output, start.elapsed());
        }

        ClassificationCommands::Isodata {
            input, output, classes, min_classes, max_classes, max_iter,
        } => {
            let raster = read_dem(&input)?;
            let start = Instant::now();
            let pb = spinner("ISODATA clustering...");
            let result = isodata(
                &raster,
                IsodataParams {
                    initial_k: classes,
                    min_k: min_classes,
                    max_k: max_classes,
                    max_iterations: max_iter,
                    min_samples: 10,
                    max_std_dev: 10.0,
                    min_merge_distance: 5.0,
                    convergence: 0.001,
                },
            )
            .context("ISODATA failed")?;
            pb.finish_and_clear();
            write_result(&result, &output, compress)?;
            done("ISODATA", &output, start.elapsed());
        }

        ClassificationCommands::Pca {
            bands, output, components,
        } => {
            let paths: Vec<&str> = bands.split(',').map(|s| s.trim()).collect();
            let pb = spinner(&format!("Reading {} bands...", paths.len()));
            let rasters: Vec<surtgis_core::Raster<f64>> = paths
                .iter()
                .map(|p| read_geotiff(p, None).with_context(|| format!("Failed to read {}", p)))
                .collect::<Result<Vec<_>>>()?;
            pb.finish_and_clear();

            let refs: Vec<&surtgis_core::Raster<f64>> = rasters.iter().collect();
            let start = Instant::now();
            let pb = spinner("PCA...");
            let result = pca(&refs, PcaParams { n_components: components })
                .context("PCA failed")?;
            pb.finish_and_clear();

            // Write each component
            let stem = output.file_stem().unwrap_or_default().to_string_lossy();
            for (i, comp) in result.components.iter().enumerate() {
                let path = output.with_file_name(format!("{}_PC{}.tif", stem, i + 1));
                write_result(comp, &path, compress)?;
                println!(
                    "  PC{}: variance={:.4} ({:.1}%) → {}",
                    i + 1,
                    result.eigenvalues[i],
                    result.variance_explained[i] * 100.0,
                    path.display()
                );
            }
            done("PCA", &output, start.elapsed());
        }

        ClassificationCommands::MinDistance {
            input, output, training,
        } => {
            let raster = read_dem(&input)?;
            let training_raster = read_dem(&training)?;
            let start = Instant::now();
            let pb = spinner("Computing signatures...");
            let sigs = signatures_from_training(&training_raster, &raster)
                .context("Failed to extract signatures")?;
            pb.finish_and_clear();
            println!("  {} class signatures extracted", sigs.len());

            let pb = spinner("Minimum distance classification...");
            let result = minimum_distance(&raster, &sigs)
                .context("Minimum distance failed")?;
            pb.finish_and_clear();
            write_result(&result, &output, compress)?;
            done("Minimum Distance", &output, start.elapsed());
        }

        ClassificationCommands::MaxLikelihood {
            input, output, training,
        } => {
            let raster = read_dem(&input)?;
            let training_raster = read_dem(&training)?;
            let start = Instant::now();
            let pb = spinner("Computing signatures...");
            let sigs = signatures_from_training(&training_raster, &raster)
                .context("Failed to extract signatures")?;
            pb.finish_and_clear();
            println!("  {} class signatures extracted", sigs.len());

            let pb = spinner("Maximum likelihood classification...");
            let result = maximum_likelihood(&raster, &sigs)
                .context("Maximum likelihood failed")?;
            pb.finish_and_clear();
            write_result(&result, &output, compress)?;
            done("Maximum Likelihood", &output, start.elapsed());
        }
    }
    Ok(())
}