Unsupervised Classification
Theory & Current API Implementation

1. Current Feature Scope

This document describes the implemented behavior, not a proposed redesign. The current registration uses analysis ID raster.classification_unsupervised, takes one multiband raster as input, allows optional band selection, and outputs one integer class raster.

2. What Unsupervised Classification Is

Unsupervised classification groups pixels using similarity in predictor-feature space without requiring labelled training polygons/classes. The algorithm discovers statistical or topological groups first. A human analyst can later interpret those groups as land-cover or material classes.

3. Current API Processing Flow

The current worker logic follows this sequence:

1. Resolve the raster input into the analysis workspace.
2. Open the raster with Rasterio and validate it.
3. Select requested predictor bands; if the band list is empty, use all bands.
4. Read selected bands as Float32.
5. Create a validity mask from finite values and raster NoData.
6. Flatten valid pixels from raster shape [bands, rows, cols] into feature matrix [pixels, bands].
7. Randomly sample at most the configured training sample size.
8. Optionally standardize predictor bands using statistics fitted on the training sample.
9. Fit the selected unsupervised algorithm.
10. Predict all valid pixels in configurable batches.
11. Add +1 to model cluster IDs so valid classes start at 1.
12. Keep invalid pixels as the configured output NoData value, default 0.
13. Write one COG with band description class.

4. Pixel Feature Space

The algorithms cluster pixels based on the selected raster values. Spatial neighbourhood is not directly used by these six algorithms in the current tool. Two distant pixels with similar selected-band vectors can therefore receive the same class.

4.1 Euclidean feature distance

K-Means, current ISODATA centre assignment, SOM BMU search, and Fuzzy C-Means distance calculations are based on Euclidean geometry.

Equivalent squared distance is often used for efficiency:

4.2 Why band selection matters

If selected bands contain redundant, noisy, saturated, cloud-contaminated, or very differently scaled variables, the clusters can be dominated by those dimensions. The tool therefore supports explicit predictor-band selection.

5. Predictor Standardization

The current default is Standardize Predictor Bands = true. The scaler is fitted on the training sample and then reused for prediction of the complete valid raster.

5.1 Why this changes clustering

Suppose two predictors have these ranges:

Band A: 0.00 – 1.00
Band B: 0 – 10,000

Without scaling, Euclidean distance is numerically dominated by Band B. After standardization, each band is approximately centered at zero with unit standard deviation, so algorithm thresholds have more consistent meaning across predictors.

6. Training Pixel Sampling & Full-Raster Prediction

The current implementation does not fit clustering on every valid raster pixel by default. It chooses a random sample, fits the model on that sample, then classifies the complete valid raster.

6.1 Fit sample

Current default:

Training Pixel Sample = 100,000
Random State = 42

6.2 Random sampling

Sampling is without replacement:

6.3 Prediction batches

The complete valid feature matrix is predicted in chunks:

Prediction Batch Size = 250,000 pixels (default)
Maximum exposed value = 2,000,000

7. MiniBatch K-Means

The current K-Means option uses MiniBatchKMeans, not the full classic K-Means implementation. The optimization target is still the within-cluster sum of squared Euclidean distances, but parameter updates are performed using mini-batches for scalability.

7.1 Objective

7.2 Assignment step

7.3 Classical center update concept

MiniBatch K-Means approximates repeated full-data updates using subsets of observations. This reduces memory and compute requirements for large pixel samples.

7.4 Current parameters

7.5 Current metadata

The analysis stores final class count and cluster centers. If standardization was enabled, centers are inverse-transformed back to the original predictor units before being stored in analysis metadata.

8. ISODATA

The current iTSensing ISODATA is a custom adaptive clustering implementation. Unlike fixed-K K-Means, it can remove small clusters, merge close centers, and split high-variance clusters. Therefore the final number of output classes may differ from the initially requested number.

8.1 Initialization

Initial centers are generated with standard K-Means, constrained between minimum and maximum classes. The internal initialization uses five K-Means starts and up to 50 iterations.

8.2 Nearest-center assignment

8.3 Center recomputation

8.4 Small-cluster removal

Clusters with fewer than min_cluster_samples are removed, unless doing so would reduce the model below the configured minimum class count. If the threshold is too aggressive, the implementation preserves the largest clusters.

8.5 Merge rule

The current implementation repeatedly finds the closest pair of cluster centers:

If:

the two centers are replaced by their simple midpoint:

8.6 Split rule

For each cluster, the implementation measures standard deviation by predictor band and chooses the band with the maximum spread:

A cluster is eligible for splitting when:

• its maximum standard deviation exceeds the split threshold;
• it has at least max(2 × min_cluster_samples, 4) members;
• splitting does not exceed the maximum class count.

The new centers are placed on the highest-variance axis:

8.7 Convergence

When the number of centers does not change, center movement is measured as:

The implementation converges only when:

8.8 Current parameters

9. Self-Organizing Map (SOM)

The current SOM learns a 2-D grid of spectral codebook vectors. After SOM training, the codebook units are grouped into the requested number of final classes using K-Means. Raster pixels are mapped to their Best Matching Unit (BMU), then inherit the class assigned to that SOM unit.

9.1 SOM map size

Current default: 5 × 5 = 25 codebook units.

9.2 Best Matching Unit

9.3 Gaussian neighbourhood

9.4 Current mini-batch target update

For each mini-batch, the neighbourhood-weighted target for a SOM unit is conceptually:

The weight update is:

9.5 Current learning-rate decay

9.6 Current neighbourhood decay

where:

9.7 Final codebook grouping

9.8 Current parameters

10. Gaussian Mixture Model (GMM)

GMM models the spectral population as a weighted mixture of K multivariate Gaussian distributions. Unlike K-Means, each component has an explicit covariance structure and therefore can represent elongated or anisotropic clusters.

10.1 Mixture density

with:

10.2 Multivariate Gaussian

10.3 Posterior responsibility

GMM fitting uses Expectation-Maximization (EM) in the underlying implementation. The current raster output uses predict(), meaning each valid pixel is assigned to one component.

10.4 Covariance types exposed by the current UI

10.5 Covariance regularization

10.6 Current parameters

10.7 Current fit metadata

The wrapper records:

• converged;
• iterations;
• lower bound;
• number of components;
• covariance type;
• cluster centers / Gaussian means.

11. Fuzzy C-Means

Fuzzy C-Means allows each training sample to have membership in every cluster. A pixel can be, for example, 0.65 related to Cluster A and 0.35 related to Cluster B during the fuzzy calculation. However, the current final raster is crisp: the implementation chooses the cluster with maximum membership.

11.1 Euclidean distance

11.2 Membership equation used by current code

Define:

Then:

where m is the fuzziness parameter and must be greater than 1.

11.3 Center update

11.4 Objective

11.5 Convergence

Stop when:

11.6 Current hard prediction

11.7 Current parameters

12. BIRCH

BIRCH (Balanced Iterative Reducing and Clustering using Hierarchies) is designed to summarize large datasets into compact Clustering Feature (CF) subclusters and then perform global clustering over those summaries. The current implementation wraps scikit-learn BIRCH and is fitted on the sampled raster feature matrix.

12.1 Clustering Feature concept

A CF summary for n observations is commonly represented by:

12.2 Centroid from CF

12.3 Current threshold

The exposed BIRCH Threshold controls the radius boundary used while building clustering features/subclusters. A lower threshold generally creates more, tighter subclusters; a higher threshold permits broader summaries.

12.4 Branching factor

The branching factor limits how many CF entries/nodes the tree can maintain before restructuring, thereby affecting memory and tree structure.

12.5 Current parameters

12.6 Current metadata

The wrapper records:

• requested clusters;
• number of generated subclusters;
• threshold;
• branching factor;
• subcluster centers.

13. Exact Current Analysis Parameters

13.1 Common parameters

13.2 Algorithm-specific visibility

The analysis registration uses dependency metadata so algorithm-specific controls appear only when the corresponding method is selected. This keeps one analysis ID while supporting six substantially different models.

14. Output Raster Contract

14.1 Raster shape

Input:
B bands × Height × Width

Output:
1 band × Height × Width

14.2 Data type

The current working output array is Int32.

14.3 NoData

Default output NoData:

0 = NoData / invalid source pixel

14.4 Valid classes

The underlying model normally predicts zero-based IDs. The current analysis explicitly adds one:

Therefore:

0 = NoData
1 = Cluster 1
2 = Cluster 2
...
K = Cluster K

14.5 COG

The result is written through the shared write_cog() utility with band description:

descriptions = ["class"]

15. Analysis Metadata Produced by the Current Tool

Common metadata fields include:

{
  "algorithm": "...",
  "bands": [...],
  "requested_classes": 5,
  "training_pixels": 100000,
  "standardized": true,
  "final_classes": ...
}

15.1 Algorithm-specific metadata

16. Algorithm Comparison

17. Which Algorithm Should a User Choose?

18. Remote-Sensing Guidance

18.1 Input should be comparable across bands

For optical imagery, Surface Reflectance is generally easier to interpret than raw DN when scenes/sensors require radiometric consistency. The clustering tool itself does not perform atmospheric correction.

18.2 Clouds and shadows

Unmasked clouds, cirrus, cloud shadows, snow, glint, and NoData artifacts can become their own clusters or distort center positions. Mask them before clustering when they are not the intended target.

18.3 Indices as predictors

The current tool expects raster bands from one input raster. If a workflow wants NDVI, NDWI, texture, DEM or other derived variables together, they should first be combined into a predictor stack using the platform's stacking/model-builder workflow.

18.4 Hyperspectral data

Hundreds of strongly correlated bands can make clustering expensive and unstable. For hyperspectral imagery, a workflow such as bad-band removal → smoothing → MNF/PCA → unsupervised clustering can be more efficient, depending on the scientific goal.

18.5 Spatial salt-and-pepper

These methods cluster spectral vectors independently; they do not enforce local spatial smoothness. A thematically noisy map is therefore possible even if the spectral clustering is mathematically valid. Post-classification majority filtering or object-based analysis is a separate operation.

19. Current Limitations

• All six methods operate on pixel feature vectors, not explicit spatial neighbourhoods.
• The model is fitted on a random sample rather than all pixels.
• There is no automatic class naming/semantic interpretation.
• There is no built-in automatic optimal-K selection such as silhouette, Calinski-Harabasz, Davies-Bouldin, BIC/AIC selection, or elbow search in this registered tool.
• Current Fuzzy C-Means final output does not expose membership bands.
• Current GMM final output does not expose posterior probability bands.
• SOM final output is a K-Means grouping of trained SOM units, not the raw 2-D SOM unit ID map.
• ISODATA split/merge thresholds depend strongly on scaling; standardization is therefore important.
• BIRCH exposes subcluster centers in metadata, but the final raster uses its requested global class prediction.
• Class IDs are arbitrary cluster identifiers and can change when parameters/random state change.

20. QA / Validation Checklist

20.1 Input validation already represented in current code

• Raster must contain valid pixels.
• Valid pixel count must be at least the requested class count.
• Training sample size must be at least the requested class count.
• Algorithm name must be one of the six supported values.
• ISODATA initial classes must fall inside minimum/maximum bounds.
• SOM must contain at least as many map units as output classes.
• SOM must have enough training samples for all units.
• GMM/Fuzzy/BIRCH require enough samples for requested components/classes.

20.2 Recommended analyst QC after classification

• Inspect number of pixels per cluster.
• Inspect cluster-center spectra/statistics.
• Compare cluster map with true-colour/false-colour imagery.
• Check whether clouds/shadows dominate clusters.
• Repeat with different band sets and class counts.
• For GMM, verify convergence.
• For ISODATA, inspect initial vs final cluster count and convergence.
• For Fuzzy C-Means, inspect objective/convergence metadata.
• Interpret/merge clusters into semantic thematic classes only after review.

21. Example Workflows

21.1 Basic Sentinel-2 land-cover exploration

Sentinel-2 Surface Reflectance
      ↓
Cloud mask
      ↓
Select B2, B3, B4, B8, B11, B12
      ↓
Unsupervised Classification
Algorithm = K-Means
Classes = 8
Standardize = true
      ↓
8-cluster COG
      ↓
Interpret / merge clusters into land-cover classes

21.2 Adaptive ISODATA

Input multispectral raster
      ↓
Classes = 8 initial
Min Classes = 4
Max Classes = 14
Min Pixels/Cluster = 50
Merge Threshold = 0.5
Split Std Threshold = 1.0
      ↓
ISODATA
      ↓
Final classes may be 6, 9, 11, etc.
      ↓
1-band class raster

21.3 GMM for overlapping spectral classes

Surface Reflectance
      ↓
Standardize = true
      ↓
GMM
Components = 6
Covariance = full
      ↓
Gaussian-mixture fit
      ↓
Maximum-posterior component ID
      ↓
6-class raster

21.4 Fuzzy C-Means

Reflectance / feature stack
      ↓
Fuzzy C-Means
m = 2.0
Classes = 5
      ↓
Soft memberships during model fitting
      ↓
argmax membership
      ↓
1-band hard 5-class raster

21.5 Hyperspectral exploratory clustering

Hyperspectral Surface Reflectance
      ↓
Bad-band removal
      ↓
Spectral smoothing
      ↓
MNF → 15 informative components
      ↓
Unsupervised Classification
K-Means / GMM / ISODATA
      ↓
Spectral group map

22. Source-to-Feature Map

This documentation is aligned to the supplied/current source modules: