Remote sensing classification and risk modeling

Updated 22 Oct 2026

•
Unsupervised classification and risk modeling
- •
  K-means clustering in remote sensing
  - •
    Workflow
    
    •
    
    •
    feature selection: select calibrated surface reflectance bands (e.g. Landsat B1-B7) and spatial region of interest.
    
    •
    sampling: extract representative training pixels across the scene.
    
    •
    centroid fitting: fit $K$ cluster centers minimizing within-cluster variance $\sum_i \|x_i - \mu_{c_i}\|^2$ .
    
    •
    pixel assignment: assign every raster pixel to its nearest cluster centroid in feature space.
    
    •
    post-classification labeling: interpret clusters into land cover classes using spectral signatures.
    
    •
    validation: evaluate accuracy against ground truth or high-resolution imagery.
  - •
    Objecive and sensitivity
    
    •
    minimizes squared Euclidean distance in feature space: $\sum_i \|x_i - \mu_{c_i}\|^2$ .
    
    •
    depends heavily on $K$ , feature scaling, initialization, and spatial sampling.
    
    •
    choosing $K$ : small $K$ merges distinct land cover classes; large $K$ splits single classes into multiple spectral subclasses.
  - •
    Interpretation cautions
    
    •
    spectral clusters are mathematical groupings, not inherent semantic categories.
    
    •
    built-up areas, roofs, and bare soil often exhibit overlapping spectral signatures.
    
    •
    always inspect all cluster centers and validate against independent reference data.
- •
  Multitemporal change detection
  - •
    Archive analysis -> track long-term land cover transitions across multi-decade satellite archives (e.g. Landsat 4/5/7 to Landsat 8/9).
  - •
    Harmonized processing workflow
    
    •
    harmonize sensor band designations across satellite generations.
    
    •
    apply cloud and shadow QA bitmasks to exclude compromised observations.
    
    •
    generate annual composites using temporal reducers (e.g. median).
    
    •
    apply calibrated spectral index classification rules (e.g. $\text{NDBI} > -0.1 \land \text{MNDWI} \le 0$ ).
  - •
    Earliest-detection raster mapping
    
    •
    assign year values to target pixels and mask non-target pixels.
    
    •
    apply a temporal minimum reducer (.reduce(ee.Reducer.min())) across the time series to identify the earliest year a pixel was observed as built-up.
    
    •
    semantic caution: represents the earliest recorded satellite detection in the series, not proof of construction date or continuous occupancy.
  - •
    Metric physical area estimation
    
    •
    Area formula -> $\text{Area (ha)} = \frac{1}{10{,}000}\sum_i \mathbf{1}_{\mathrm{built}, i} \cdot A_i$ , where $A_i$ is metric pixel area in $\text{m}^2$ .
    
    •
    Calculation rule -> sum binary target indicator masks; do not multiply area by raster category codes or year values.
  - •
    Temporal comparability cautions
    
    •
    inter-annual variations in cloud cover, seasonal timing, and sensor characteristics can cause artificial area fluctuations.
    
    •
    date intervals in Earth Engine: filterDate with an end date of December 31 is end-exclusive; use January 1 of the following year to include all December observations.

Python implementation: clustering and multitemporal change detection

•
demonstrates local k-means clustering, earliest-detection temporal aggregation, and correct metri area computation. uses require_plane from code 3.

import cv2
import numpy as np

def cluster_features(features, k=2, sample_size=500, seed=0):
        a = np.asarray(features, dtype=np.float32)
        if a.ndim != 3 or a.size == 0 or not np.isfinite(a).all():
            raise ValueError("expected finite H x W x bands features, already validity-filtered")
        pixels = a.reshape(-1, a.shape[-1])
        if not isinstance(k, int) or not isinstance(sample_size, int) or not 1 <= k <= sample_size:
            raise ValueError("require integer 1 <= k <= sample_size")
        if k > len(pixels) or not isinstance(seed, int) or not 0 <= seed < 2**31:
            raise ValueError("not enough pixels or unsupported seed")
        rng = np.random.default_rng(seed)
        selected = rng.choice(len(pixels), min(sample_size, len(pixels)), replace=False)
        if len(selected) == 1:
            # OpenCV treats a single-row matrix as scalar samples, not one feature vector.
            centers = pixels[selected].copy()
        else:
            cv2.setRNGSeed(seed)  # OpenCV global RNG: record seed; avoid concurrent reseeding
            criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_MAX_ITER, 100, 1e-4)
            _, _, centers = cv2.kmeans(pixels[selected], k, None, criteria, 1, cv2.KMEANS_PP_CENTERS)
        distances = ((pixels[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
        return distances.argmin(axis=1).reshape(a.shape[:2]), centers

def temporal_built_summary(built, valid, years, pixel_area_m2):
        if not isinstance(built, np.ndarray) or built.dtype != bool or built.ndim != 3 or built.size == 0:
            raise ValueError("built must be nonempty Boolean T x H x W")
        if not isinstance(valid, np.ndarray) or valid.dtype != bool or valid.shape != built.shape:
            raise ValueError("valid must match built")
        years = np.asarray(years)
        if years.ndim != 1 or len(years) != built.shape[0] or not np.issubdtype(years.dtype, np.integer):
            raise ValueError("one integer year per time step required")
        if np.any(years <= 0) or np.any(np.diff(years) <= 0):
            raise ValueError("years must be positive and strictly increasing")
        area = require_plane(pixel_area_m2)
        if area.shape != built.shape[1:] or np.any(area <= 0):
            raise ValueError("positive aligned pixel areas required")
        detected = built & valid
        first = np.min(np.where(detected, years[:, None, None], np.inf), axis=0)
        first[~detected.any(axis=0)] = np.nan  # never detected, not year zero
        annual_ha = (detected * area[None]).sum(axis=(1, 2)) / 10000
        coverage_ha = (valid * area[None]).sum(axis=(1, 2)) / 10000
        return first, annual_ha, coverage_ha

features = np.array([[[0., 0.], [0., 0.]], [[1., 1.], [1., 1.]]], np.float32)
labels, centers = cluster_features(features, k=2, sample_size=4, seed=4)
assert labels[0, 0] == labels[0, 1] and labels[0, 0] != labels[1, 0]
single_labels, single_center = cluster_features(
        np.array([[[0., 10.]]]), k=1, sample_size=1)
assert single_center.shape == (1, 2) and np.array_equal(single_center, [[0., 10.]])
assert single_labels.tolist() == [[0]]
built = np.array([[[True, False]], [[False, True]]])
first, annual_ha, coverage_ha = temporal_built_summary(
        built, np.ones_like(built), [1990, 2000], np.full((1, 2), 900.))
assert first.tolist() == [[1990., 2000.]]
assert np.allclose(annual_ha, [0.09, 0.09]) and np.allclose(coverage_ha, [0.18, 0.18])

•
Multi-criteria environmental risk modeling
- •
  Hazard, Exposure, and Vulnerability (HEV) framework
  - •
    Risk definition -> disaster risk is modeled as the interaction of hazard, exposure, and vulnerability ( $R = f(H, E, V)$ ).
  - •
    Hazard ( $H$ ) -> physical likelihood, magnitude, and spatial extent of dangerous events (e.g. flood inundation).
  - •
    Exposure ( $E$ ) -> population, assets, infrastructure, and land cover located in hazard zones.
  - •
    Vulnerability ( $V$ ) -> socioeconomic susceptibility to harm and capacity to anticipate, cope, and recover.
- •
  Harkat scoring (ordinal classification)
  - •
    reclassifies heterogeneous continuous variables into a common ordinal scale (typically scores 1 to 5).
  - •
    Ascending score -> higher values correspond to higher classes: $s_\uparrow(x) = 1 + \sum_j \mathbf{1}(x > b_j)$ .
  - •
    Descending score -> lower values correspond to higher hazard/risk: $s_\downarrow(x) = 6 - \left[1 + \sum_j \mathbf{1}(x > b_j)\right]$ .
- •
  Factor selection and derivation
  - •
    Hazard components ( $H$ )
    
    •
    distance to water bodies: proximity to river channels increases flood potential (descending score).
    
    •
    elevation: low elevation accumulates runoff (descending score).
    
    •
    Topographic Position Index (TPI): relative topographic depressions (negative TPI) pool water.
    
    •
    Vegetation cover (NDVI): low vegetation or high impervious cover accelerates surface runoff.
    
    •
    terrain slope: flat terrain impedes drainage, heightening pooling hazard.
  - •
    Exposure components ( $E$ )
    
    •
    population density grids (e.g. GHSL population counts per cell).
    
    •
    land cover classification: built-up human settlements and agricultural crop fields.
  - •
    Vulnerability components ( $V$ )
    
    •
    socioeconomic capacity proxies: GDP per capita and Human Development Index (HDI).
    
    •
    lower economic resources and lower HDI reduce capacity to absorb and recover from losses.
- •
  Risk index aggregation
  - •
    Additive combination
    
    •
    two-component model: $r_2 = H + E$ .
    
    •
    three-component model: $r_3 = H + E + V$ .
  - •
    Weighted multi-criteria model -> $R = w_H H + w_E E + w_V V$ , where weights $\sum w_i = 1$ reflect relative influence.
  - •
    Scale comparability
    
    •
    comparing raw sums across models with different numbers of components ( $r_2 \in [2, 10]$ vs $r_3 \in [3, 15]$ ) introduces range artifacts.
    
    •
    normalize by component count ( $r_2/2$ vs $r_3/3$ ) or adopt standardized scale cuts (e.g. 1.8, 2.6, 3.4, 4.2) for meaningful comparative analysis.
- •
  Sensitivity analysis and model limitations
  - •
    Sensitivity testing -> evaluate model stability against variations in factor weights and classification cutoffs.
  - •
    Index vs hydraulic modeling -> multi-criteria harkat scoring provides relative spatial rankings of indicative risk; it does not simulate fluid mechanics or output calibrated physical flood depths/probabilities.
  - •
    Empirical validation -> validate predicted risk hotspots against historical flood inundation footprints and damage records.

•

Python implementation: multi-criteria risk modeling

•
computes ordinal harkat scores, combines multi-criteria risk models, and summarizes class areas in hectares. reuses require_plane form code 3

•

import numpy as np

def harkat(values, cuts, descending=False):
        x = require_plane(values)
        cuts = np.asarray(cuts, dtype=np.float64)
        if cuts.shape != (4,) or not np.isfinite(cuts).all() or np.any(np.diff(cuts) <= 0):
            raise ValueError("four finite strictly increasing cuts required")
        score = 1 + (x[..., None] > cuts).sum(axis=-1)
        return (6 - score if descending else score).astype(np.uint8)

def flood_model(distance, dem, tpi, ndvi, exposure, vulnerability, slope, valid,
                    slope_weight=0.2, slope_cuts=(2, 8, 15, 30)):
        planes = list(map(require_plane, (distance, dem, tpi, ndvi, exposure, vulnerability, slope)))
        if len({x.shape for x in planes}) != 1:
            raise ValueError("all model variables must have the same grid")
        distance, dem, tpi, ndvi, exposure, vulnerability, slope = planes
        if not isinstance(valid, np.ndarray) or valid.dtype != bool or valid.shape != dem.shape:
            raise ValueError("a common Boolean validity mask is required")
        if np.any(distance < 0) or np.any(slope < 0):
            raise ValueError("distance and slope cannot be negative")
        if any(np.any((x < 1) | (x > 5) | (x != np.rint(x))) for x in (exposure, vulnerability)):
            raise ValueError("E and V must be previously classified integer scores 1..5")
        if not np.isfinite(slope_weight) or not 0 <= slope_weight <= 1:
            raise ValueError("slope weight must be in 0..1")
        scores = [harkat(x, cuts, descending=True) for x, cuts in (
            (distance, (50, 100, 150, 200)), (dem, (5, 10, 15, 20)),
            (tpi, (-0.75, -0.5, -0.25, 0)), (ndvi, (-0.25, 0, 0.25, 0.5)))]
        hazard_sum = np.stack(scores).astype(np.float64).sum(axis=0)
        H = harkat(hazard_sum, (4, 8, 12, 16))
        r2, r3 = H + exposure, H + exposure + vulnerability
        slope_score = harkat(slope, slope_cuts, descending=True)
        new_hazard = hazard_sum / 4 * (1 - slope_weight) + slope_score * slope_weight
        weighted = 0.5 * new_hazard + 0.3 * exposure + 0.2 * vulnerability
        outputs = {"H": H, "r2": r2, "r3": r3, "weighted": weighted,
                   "k2": harkat(r2, (2, 4, 6, 8)),
                   "k3_literal": harkat(r3, (2, 4, 6, 8)),
                   "k3_table": harkat(r3, (2, 6, 9, 12)),
                   "k2_mean": harkat(r2 / 2, (1.8, 2.6, 3.4, 4.2)),
                   "k3_mean": harkat(r3 / 3, (1.8, 2.6, 3.4, 4.2)),
                   "weighted_class": harkat(weighted, (1.8, 2.6, 3.4, 4.2))}
        return {name: np.where(valid, x, np.nan) for name, x in outputs.items()}

def class_area_ha(classes, pixel_area_m2):
        classes = np.asarray(classes, dtype=np.float64)
        area = require_plane(pixel_area_m2)
        if classes.shape != area.shape or np.isinf(classes).any() or np.any(area <= 0):
            raise ValueError("aligned classes and positive pixel areas required")
        ok = np.isfinite(classes)
        if np.any((classes[ok] < 1) | (classes[ok] > 5) | (classes[ok] != np.rint(classes[ok]))):
            raise ValueError("valid classes must be integer labels 1..5; NaN means unavailable")
        return {k: float(area[classes == k].sum() / 10000) for k in range(1, 6)}

cuts_example = np.array([[50., 100., 150., 200., 201.]])
assert harkat(cuts_example, (50, 100, 150, 200), True).tolist() == [[5, 4, 3, 2, 1]]
plane = np.ones((1, 2))
risk = flood_model(plane * 250, plane * 25, plane, plane * 0.6,
                       plane, plane, plane * 35, np.array([[True, False]]))
assert risk["r2"][0, 0] == 2 and risk["r3"][0, 0] == 3
assert risk["k3_literal"][0, 0] == 2 and risk["k3_mean"][0, 0] == 1
assert risk["weighted"][0, 0] == 1 and np.isnan(risk["weighted"][0, 1])
assert class_area_ha(risk["weighted_class"], plane * 900)[1] == 0.09

•
additive ordinal scores and selected weights form an indicative ranking
- •
  does not provide calibrated flood probability, hydraulic simulation, or guarantee of flooding
- •
  validate against actual events
- •
  rainfall, drainage, and source uncertainty are other infuences

Linked references 1

digital-image-processing-polban

◌ Explore connections in Graph view

Remote sensing classification and risk modeling

Unsupervised classification and risk modeling

Python implementation: clustering and multitemporal change detection

Multi-criteria environmental risk modeling

Python implementation: multi-criteria risk modeling

Linked references 1