🧭 Concept Natural hazards & risk Terrain & elevation Natural Hazards & Risk · step 6 Notebook By · Updated

Landslide susceptibility explained: what the map is and is not

On this page

Problem statement

A landslide susceptibility map ranks the ground by how prone it is to failure. It does not say when, it does not say how often, and unless it was built and tested against an inventory of landslides that actually happened, it does not say how well it works either. Those three absences are what separate susceptibility from hazard and from risk — and they are usually missing from the legend.

The ingredients are easy to assemble: slope from a DEM, land cover, distance to streams, maybe lithology. Combining them into a number is easy too. What is hard, and what decides whether the map is worth anything, is an inventory of past failures to fit the weights and to measure the result against.

Quick answer

Build the factors, combine them, and be explicit about which kind of map you have made:

import numpy as np, rasterio, pystac_client, planetary_computer
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
from sar_io import STAC, grid                      # the grid helper from the radar cluster

bbox = [-122.00, 40.00, -121.74, 40.22]            # steep ground inside the 2024 Park Fire
g = grid(bbox, "EPSG:32610", res=30)
cat = pystac_client.Client.open(STAC, modifier=planetary_computer.sign_inplace)

def on_grid(item, asset, resampling=Resampling.bilinear):
    with rasterio.open(item.assets[asset].href) as src, WarpedVRT(
            src, crs=g["crs"], transform=g["transform"], width=g["width"],
            height=g["height"], resampling=resampling, nodata=np.nan) as vrt:
        return vrt.read(1)

dem = np.full((g["height"], g["width"]), np.nan, "float32")
for it in cat.search(collections=["cop-dem-glo-30"], bbox=bbox).items():
    a = on_grid(it, "data")
    dem[np.isnan(dem) & np.isfinite(a)] = a[np.isnan(dem) & np.isfinite(a)]

gy, gx = np.gradient(dem, 30)
slope = np.degrees(np.arctan(np.hypot(gx, gy)))
print(f"slope: median {np.nanmedian(slope):.1f}°, 95th percentile {np.nanpercentile(slope, 95):.1f}°, "
      f"steeper than 30°: {np.nanmean(slope > 30):.1%} of the area")
slope: median 13.9°, 95th percentile 25.4°, steeper than 30°: 2.2% of the area

Slope is the one factor that appears in every susceptibility model ever published. It is also, on its own, not a model.

Susceptibility, hazard and risk for landslides: what each adds and what each can answer. three products, three questionsproductaddscan answersusceptibilitya ranking of placeswhere is it steeper trouble?hazard+ frequency and sizehow often, how big?risk+ exposure, vulnerabilitywhat will it cost?most published maps are the first row; most decisions need the thirdA susceptibility class is relative to the rest of the map, and to nothing else spatialworkflow.io
Three products, three questions; only the last one supports a decision about a building.

Physical factors, an inventory, honest absences, spatial validation

1. Choose factors with a physical reason

Slope, slope curvature, aspect, land cover, distance to streams and roads, lithology, soil depth. Each should have a mechanism behind it, not just a correlation in one study area.

2. Get an inventory

Mapped landslide polygons or points, ideally with dates. Without an inventory there is nothing to fit weights to and nothing to test against.

3. Choose a method

Heuristic weighting by an expert; statistical fitting — weights of evidence, logistic regression, random forest — against the inventory; or physically based slope-stability modelling where soil and hydrology data exist.

4. Sample absences honestly

Statistical models need non-landslide locations as well as landslide ones. Where those come from decides the result: sampling them from flat valley floors guarantees a model that has learned "steep is bad" and nothing else.

5. Validate spatially

Hold out an area, not a random subset of pixels, and report the area under the ROC curve on the held-out ground. Landslide inventories are strongly clustered, so a random split is far too easy.

6. Label the product correctly

Susceptibility is relative propensity with no time attached. Hazard adds a frequency. Risk adds exposure and vulnerability.

Share of a study area in each susceptibility class of a quantile-classified index, with median slope. quantile classes of a heuristic indexmoderate or above49.7% · median slope 18°high or above19.9% · 22°very high (top 5%)5.0% · 29°steeper than 30°2.2% of the area050%100%Quantile classes cover a fixed share of the map by construction, whatever the ground is like spatialworkflow.io
A susceptibility class covers a share of the map; without an inventory, that share is all you know about it.

A heuristic index and what validating it would need

The data used on this page

The terrain is the Copernicus GLO-30 DEM from Microsoft Planetary Computer, read onto a 30 m UTM grid over steep ground inside the 2024 Park Fire in northern California. sar_io.py is the small helper module from the radar guides; this page uses only its STAC address and its grid() function.

%%writefile sar_io.py
"""Find, read and mosaic Sentinel-1 RTC scenes from Microsoft Planetary Computer."""
from collections import defaultdict

import numpy as np
import planetary_computer
import pystac_client
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
from rasterio.warp import transform_bounds

STAC = "https://planetarycomputer.microsoft.com/api/stac/v1"


def search(bbox, dates, orbit=None, direction=None, collection="sentinel-1-rtc"):
    """Items over bbox (lon/lat) in a date range, filtered by relative orbit and pass direction."""
    cat = pystac_client.Client.open(STAC, modifier=planetary_computer.sign_inplace)
    items = cat.search(collections=[collection], bbox=bbox, datetime=dates).item_collection()
    keep = [i for i in items
            if (orbit is None or i.properties["sat:relative_orbit"] == orbit)
            and (direction is None or i.properties["sat:orbit_state"] == direction)]
    return sorted(keep, key=lambda i: i.datetime)


def passes(items):
    """Group frames into passes: same day, same relative orbit."""
    groups = defaultdict(list)
    for i in items:
        groups[(i.datetime.date(), i.properties["sat:relative_orbit"])].append(i)
    return dict(sorted(groups.items()))


def grid(bbox, crs, res=10):
    """A fixed target grid: every pass is read onto this, so arrays line up."""
    left, bottom, right, top = transform_bounds(4326, crs, *bbox)
    width, height = int((right - left) // res), int((top - bottom) // res)
    return {"crs": crs, "transform": from_origin(left, top, res, res),
            "width": width, "height": height}


def read_pass(frames, band, g):
    """Mosaic one pass's frames for one band onto grid g. Linear power; NaN where no data."""
    out = np.full((g["height"], g["width"]), np.nan, dtype="float32")
    for item in frames:
        with rasterio.open(item.assets[band].href) as src, WarpedVRT(
                src, crs=g["crs"], transform=g["transform"], width=g["width"],
                height=g["height"], resampling=Resampling.bilinear, nodata=np.nan) as vrt:
            a = vrt.read(1)
        a[a <= 0] = np.nan                       # -32768 nodata and zero-padding
        fill = np.isnan(out) & np.isfinite(a)
        out[fill] = a[fill]
    return out


def to_db(linear):
    return 10 * np.log10(linear)
Writing sar_io.py

Example 1 — a heuristic index, and what it can and cannot say

wc_item = next(cat.search(collections=["esa-worldcover"], bbox=bbox).items())
wc = on_grid(wc_item, "map", Resampling.nearest)

def rescale(x, lo, hi):
    return np.clip((x - lo) / (hi - lo), 0, 1)

slope_f = rescale(slope, 10, 45)                       # steeper is worse, saturating at 45°
cover_f = np.select([np.isin(wc, [50, 60]), np.isin(wc, [30, 40]), np.isin(wc, [10, 20])],
                    [0.9, 0.6, 0.3], default=0.5)      # bare and built-up worst, forest best
index = 0.7 * slope_f + 0.3 * cover_f

qs = np.nanpercentile(index, [50, 80, 95])
print("index percentiles:", np.round(qs, 2).tolist())
for name, lo in [("moderate", qs[0]), ("high", qs[1]), ("very high", qs[2])]:
    m = index >= lo
    print(f"{name:10} {np.nanmean(m):5.1%} of the area, median slope {np.nanmedian(slope[m]):4.1f}°")
index percentiles: [0.22, 0.32, 0.44]
moderate   49.7% of the area, median slope 18.0°
high       19.9% of the area, median slope 22.0°
very high   5.0% of the area, median slope 29.4°

The classes here are quantiles of the index — by construction the top class is 5% of the area. That is a statement about the map, not about the ground.

Example 2 — what validation would need

def auc(scores, labels):
    """Area under the ROC curve, computed from ranks: the probability that a
    randomly chosen landslide pixel scores above a randomly chosen stable one."""
    order = np.argsort(scores)
    ranks = np.empty(len(scores)); ranks[order] = np.arange(1, len(scores) + 1)
    pos = labels == 1
    n1, n0 = pos.sum(), (~pos).sum()
    return (ranks[pos].sum() - n1 * (n1 + 1) / 2) / (n1 * n0)

rng = np.random.default_rng(0)
fake_scores = rng.random(1000)
fake_labels = (rng.random(1000) < 0.2).astype(int)
print(f"AUC of a random score against random labels: {auc(fake_scores, fake_labels):.3f}")
print("With a real inventory, the same function takes the index values at landslide")
print("locations and at stable locations, and returns how well the map separates them.")
AUC of a random score against random labels: 0.525
With a real inventory, the same function takes the index values at landslide
locations and at stable locations, and returns how well the map separates them.

Example 3 — susceptibility is not probability

top = index >= qs[2]
print(f"area in the top class: {top.sum() * 900 / 1e6:.1f} km²")
print("annual probability of a landslide there: unknown from this map")
print("expected losses there: unknown from this map")
print("what the map does say: this ground ranks above 95% of the study area")
area in the top class: 27.5 km²
annual probability of a landslide there: unknown from this map
expected losses there: unknown from this map
what the map does say: this ground ranks above 95% of the study area

Why an inventory is not optional

Why an inventory is not optional

Weights fitted without data are opinions, and opinions cannot be tested. An inventory gives both: it fits the weights and it measures how well the resulting ranking separates places that failed from places that did not. A susceptibility map with no inventory behind it can still be useful as a screening layer, but it should say so, and it should not be given class names that imply measurement.

Why absence sampling decides the result

A statistical model needs examples of stable ground. Take them anywhere, and most will be flat, so the model learns that steep ground is dangerous — which was already known — and nothing about the rest. Sampling absences from within the same terrain envelope as the landslides forces the model to explain something harder.

Why susceptibility, hazard and risk are different words

Susceptibility ranks places. Hazard attaches a frequency: how often a failure of a given size happens there. Risk multiplies by what is exposed and how vulnerable it is. Most published maps are susceptibility, most decisions need risk, and the gap between them is usually left for the reader to notice.

Absence sampling anywhere in the study area against sampling within the same terrain envelope. where the "no landslide" examples come fromabsences sampled anywheremost land is flat✕model learns "steep is bad"✕AUC looks excellent✕adds nothing to a slope mapabsences from the same terrainsteep, stable ground included✓model must explain more✓AUC is lower and honest✓tells you what slope does notThe choice of absences decides what a statistical susceptibility model can learn spatialworkflow.io
Where the "no landslide" examples come from decides what the model can learn.

Triggers, post-fire debris flows and runout

  • Triggering matters: rainfall-triggered and earthquake-triggered landslides have different susceptibility patterns.
  • Post-fire debris flows follow burn severity as much as slope, and the window of elevated risk lasts a few years.
  • Inventory completeness varies with the imagery used to map it; small slides are missed.
  • DEM resolution changes slope: a 30 m DEM smooths away the steep faces where failures start.
  • Runout is a separate model; source susceptibility does not say where the debris stops.
  • Report the class shares; a "very high" class covering a third of the map is not informative.

FAQ

What is landslide susceptibility?

A relative ranking of how prone ground is to failure, with no time or frequency attached.

How is it different from landslide hazard?

Hazard adds a frequency — how often a failure of a given size occurs in that place.

Can I build one without an inventory?

You can build an index. You cannot fit weights or measure performance, and the map should say so.

What factors should I use?

Slope first, then curvature, land cover, distance to streams, and lithology where available.

How do I validate it?

Area under the ROC curve against a held-out area of the inventory, not a random pixel split.

Does a "very high" class mean a landslide is likely?

Only relative to the rest of the map. Without an inventory and a frequency, no probability is implied.