Short answer

Landslide susceptibility mapping estimates where slope failures are spatially more likely, by combining a complete landslide inventory with terrain, geological, and hydrological predictor layers in a statistical or heuristic model. It answers a relative question ("which slopes resemble those that have already failed?"), not an absolute one ("will this slope fail this winter?"). The two most defensible production methods are data-driven statistics — weights-of-evidence, frequency ratio, or logistic regression — validated against a withheld portion of the inventory using the Area Under the Curve (AUC). A slope raster on its own is not a susceptibility map.

Below is a concrete workflow with the predictor layers, the actual algorithms, the GDAL and QGIS steps, and the validation maths.

Susceptibility, hazard, and risk are three different products

These words are used loosely in practice, and conflating them is the most common reporting error.

  • Susceptibility is the spatial likelihood of failure given terrain and ground conditions. It has no time component and no measure of magnitude.
  • Hazard adds frequency and intensity: the probability that a landslide of a given size occurs in a given area within a given period. This requires a dated inventory and rainfall or seismic triggering data.
  • Risk multiplies hazard by exposure and vulnerability — the people, buildings, or infrastructure in the runout path and their expected loss.

A GIS susceptibility model delivers the first product. Calling it a hazard or risk map without the temporal and consequence components overstates it.

The predictor layers

Susceptibility models combine conditioning factors (slow-changing terrain and ground properties) and, for hazard, triggering factors (rainfall, seismicity). The standard conditioning set:

Terrain derivatives from a DEM

Start from a hydrologically sound DEM — a 10 m or 12.5 m DEM (for example USGS 3DEP, or ALOS PALSAR RTC) is a common compromise; 1–2 m LiDAR is better where available but can introduce noise that needs smoothing.

  • Slope angle — the dominant predictor. gdaldem slope dem.tif slope.tif -compute_edges. Most shallow translational failures cluster in a characteristic band, often roughly 25°–45°; very steep rock faces and near-flat ground both show lower shallow-failure density.
  • Aspect — proxies for solar exposure, prevailing rain direction, and bedding orientation: gdaldem aspect dem.tif aspect.tif.
  • Plan and profile curvature — concave hollows concentrate water and debris; compute with QGIS Processing gdal:slopeaspectcurvature-style tools or SAGA Slope, Aspect, Curvature.
  • Topographic Wetness Index (TWI) = ln(a / tan β), where a is the upslope contributing area per unit contour length and β is local slope. SAGA's Topographic Wetness Index algorithm computes it directly; high TWI flags saturation-prone hollows.

Geological and hydrological factors

  • Lithology / engineering geology units, reclassified to relative strength or weathering susceptibility.
  • Distance to faults / structural lineamentsgdal:proximity on a rasterised fault layer; fractured, sheared rock near faults is weaker.
  • Distance to drainage — toe erosion undercuts slopes; again use gdal:proximity.
  • Land cover / NDVI — deep-rooted forest stabilises soil; recent clearance or burn scars raise susceptibility.

The landslide inventory

This is the part teams under-invest in, and it controls everything downstream. You need polygons or points of known failures, ideally separating the source/scarp area from the deposit, because predictors should be sampled at the depletion zone, not the runout. Map from high-resolution imagery, LiDAR hillshades, or existing agency inventories. Record the failure type (shallow translational, deep-seated rotational, debris flow, rockfall) — a single model mixing all types performs poorly.

A worked statistical workflow (frequency ratio / weights-of-evidence)

  1. Harmonise the data cube. Reproject every layer to one projected CRS and resample to a common pixel size — a national grid such as EPSG:32633 (UTM) at 12.5 m is typical. gdalwarp -t_srs EPSG:32633 -tr 12.5 12.5 -r bilinear in.tif out.tif for continuous rasters; use -r near for categorical layers (lithology, land cover) so classes are not interpolated.

  2. Reclassify continuous predictors into classes. For example, slope into bins (0–10, 10–20, 20–30, 30–40, 40–60, >60). Use QGIS Reclassify by table or gdal_calc.

  3. Split the inventory. Randomly assign 70% of inventory cells to training and hold 30% for validation. Keep the split reproducible (fixed random seed).

  4. Compute frequency ratio (FR) per class. For each class: FR = (landslide pixels in class / total landslide pixels) / (total pixels in class / total pixels). FR > 1 means failures are over-represented in that class. Sum the FR values of the contributing classes per pixel to get the Landslide Susceptibility Index (LSI). Weights-of-evidence is the Bayesian formalisation of the same idea, producing positive/negative weights (W⁺, W⁻) and a contrast C that doubles as a factor-significance test.

  5. Alternatively, logistic regression. Sample predictor values at landslide and an equal number of randomly placed non-landslide points, then fit P = 1 / (1 + e^-(b0 + b1x1 + ... )). This handles correlated, continuous predictors better and yields a 0–1 probability surface. Check predictor multicollinearity first (variance inflation factor; drop or merge factors with VIF > 5–10).

  6. Classify the index into 4–5 relative classes (Very Low to Very High) using natural breaks (Jenks) or quantiles, and document which method you used — the class boundaries are an interpretive choice, not a fact.

Validating the map

Never publish a susceptibility map without a quantitative validation.

  • Overlay the withheld 30% of inventory on the susceptibility surface.
  • Build a prediction-rate curve: sort pixels from highest to lowest susceptibility, then plot cumulative percentage of validation landslides against cumulative percentage of study area.
  • Compute the AUC. As a rough guide: 0.5 = no skill, 0.7–0.8 = acceptable, 0.8–0.9 = good, >0.9 = excellent (and worth checking for inventory leakage or overfitting). The "success rate" curve (training data) measures fit; the "prediction rate" curve (test data) measures real predictive power — report both, and trust the prediction rate.

A model that puts most validation failures in the top 20–30% of ranked area is performing as intended.

Common pitfalls and why they happen

  • Incomplete or biased inventory. If older or vegetated failures are missed, the model learns the mapping bias, not the physics. Inventories skewed toward roadside failures over-weight "distance to road."
  • Sampling the runout instead of the source. Predictor values in the deposit reflect where debris stopped, not where it failed. Use scarp/source polygons.
  • Resolution mismatch. Mixing a 30 m DEM with a 1:10,000 lithology map gives a false sense of precision; the product is only as sharp as the coarsest meaningful input.
  • Spatial autocorrelation in the split. Random pixel splitting leaks information between adjacent training and test cells, inflating AUC. Split by landslide polygon or by spatial block instead.
  • Categorical resampling with bilinear interpolation. Interpolating lithology codes invents non-existent classes. Always use nearest-neighbour for class rasters.
  • Treating a relative index as failure probability. A "Very High" class means "resembles past failures," not "this slope has an 80% chance of sliding."

QA and validation checklist

  • DEM source, resolution, and date recorded; sinks filled before deriving flow-based indices.
  • All layers in one projected CRS, one pixel size, identical extent and grid alignment.
  • Inventory completeness, failure-type filtering, and source-vs-deposit handling documented.
  • Train/test split method stated (and ideally spatially blocked).
  • Predictor multicollinearity checked for regression models.
  • AUC reported for prediction rate, not just success rate.
  • Susceptibility classes labelled as relative, with the classification method named.

Bathyl perspective

We treat a susceptibility map as a screening and prioritisation tool that focuses field effort and routing decisions — not a substitute for geotechnical site investigation. The deliverable that earns trust is the one shipped with its inventory, its AUC, and an explicit statement of what failure types and time horizon it does and does not cover. When a map drives permitting, design, or public safety, it should be the entry point to specialist assessment, not the conclusion.

Related reading

Sources