Short answer

Erosion risk mapping in GIS usually means estimating long-term soil loss with the Revised Universal Soil Loss Equation (RUSLE), A = R × K × LS × C × P, where each factor is a raster aligned to a common grid and the product gives average annual soil loss in tonnes per hectare per year. The terrain part of the work, the slope-length and steepness (LS) factor, comes directly from a DEM, which is why coordinate and resolution discipline matter as much as the erosion science.

The output is a screening surface that ranks where erosion is most likely, not a guarantee of where it will occur. RUSLE captures sheet and rill erosion under average conditions; it says nothing about gully incision, streambank retreat, or landslides, and it is highly sensitive to DEM cell size. Used honestly, it tells a planner where to look first.

What RUSLE actually models (and what it ignores)

RUSLE is an empirical model fitted to plot data, mostly agricultural plots in the United States. It predicts the long-term average annual soil loss from sheet and rill erosion on a hillslope. The five factors are:

  • R — rainfall-runoff erosivity (MJ·mm·ha⁻¹·h⁻¹·yr⁻¹). Higher where storms are intense.
  • K — soil erodibility (t·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹). A property of texture, organic matter, structure, and permeability.
  • LS — slope length and steepness (dimensionless). The terrain amplifier, derived from the DEM.
  • C — cover-management (0–1). Bare soil is near 1; dense forest is near 0.
  • P — support practice (0–1). Contouring, terracing, and strip cropping lower it; 1 means no practice.

Three things RUSLE is not: it is not a sediment-yield model (it estimates gross hillslope loss, not what reaches a channel; for delivery you need a sediment delivery ratio or a model like SDR in InVEST), it does not represent concentrated-flow gully or ephemeral-gully erosion, and it does not handle mass movement. If your real concern is landslides or debris flows, RUSLE is the wrong tool and you want a susceptibility model instead.

Building the LS factor from a DEM

The LS factor is where GIS earns its place. Start with a hydrologically conditioned DEM in a projected CRS with metres (e.g. UTM, EPSG:326xx/327xx, or a national grid). Angular degrees will silently corrupt slope and length calculations.

A common workflow in QGIS or GDAL:

  1. Fill sinks so flow routing is continuous. QGIS Processing offers Fill sinks (Wang & Liu) (SAGA) or r.fill.dir (GRASS). Over-aggressive filling flattens real depressions, so inspect the difference raster.
  2. Slope in degrees: gdaldem slope dem.tif slope_deg.tif -alg Horn. The Horn algorithm is the standard 3×3 estimator.
  3. Flow accumulation: GRASS r.watershed or SAGA Flow Accumulation (Top-Down) gives the number of upslope cells, which proxies slope length.
  4. Combine into LS with the Moore & Burch (1986) formulation, a widely used GIS-friendly form:
LS = (flowacc * cellsize / 22.13)^0.4 * (sin(slope_rad) / 0.0896)^1.3

In QGIS Raster Calculator or gdal_calc.py you express slope_rad = slope_deg * pi/180. The constants 22.13 m and 0.0896 are the RUSLE unit-plot reference length and slope. Be deliberate about cell size: an LS surface from a 2 m LiDAR DEM and one from a 30 m SRTM DEM are not comparable, because flow accumulation scales with the grid.

Assembling the other factors

  • K comes from a soils polygon layer (e.g. SSURGO in the US, or a national soil map) joined to a Kfactor attribute, then rasterized to the working grid with gdal_rasterize -a kfactor -tr 10 10.
  • R is usually an interpolated isoerodent surface or a published gridded product. If you only have annual precipitation, the Renard-Freimund or Wischmeier approximations give a defensible first estimate; document which you used.
  • C can be reclassified from a land-cover raster (CORINE, ESA WorldCover, NLCD) using literature C-values, or estimated from NDVI with C = exp(-α · NDVI/(β − NDVI)). NDVI gives a continuous, season-specific cover surface but inherits any cloud or date bias from the imagery.
  • P is often set to 1 unless you have mapped conservation practices.

Reproject and resample everything to one grid first. A single gdalwarp -t_srs EPSG:32631 -tr 10 10 -te <xmin ymin xmax ymax> -r bilinear per layer (use -r near for categorical inputs like land cover) guarantees identical extent, cell size, and alignment, which is non-negotiable before map algebra.

Worked example

Suppose a 10 m grid catchment with: R = 600, a silt-loam soil K = 0.045, a cultivated parcel C = 0.35, P = 1, and a hillslope cell where LS = 4.2.

A = 600 * 0.045 * 4.2 * 0.35 * 1 = 39.7 t/ha/yr

That sits in a high-erosion class (many guidelines treat >10–20 t/ha/yr as severe for cropland). Re-run with cover restored (C = 0.04 for grassland) and A drops to about 4.5 t/ha/yr, which quantifies the value of cover and is exactly the kind of scenario comparison RUSLE supports well. The full raster version:

gdal_calc.py -R r.tif -K k.tif --L=ls.tif -C c.tif -P p.tif \
  --calc="R*K*L*C*P" --outfile=soilloss.tif --NoDataValue=-9999

Then classify soilloss.tif into erosion classes for the map, and keep the continuous raster for downstream analysis.

Common pitfalls and why they happen

  • Mixing DEM resolutions across the catchment. LS depends on flow accumulation, which is grid-dependent, so a mosaic of 5 m and 30 m tiles produces a discontinuity at the seam. Resample to one resolution first.
  • Leaving the DEM in geographic coordinates. Slope computed on degrees treats a 0.0001° cell as if it were a metre, producing nonsense LS values. This is the single most frequent error.
  • Treating RUSLE output as sediment yield. Gross soil loss overstates what reaches a stream; without a delivery ratio you will over-promise downstream impact.
  • Reusing literature C/K values without noting the source region. Factor tables are regional; copying US cropland C-values into a Mediterranean shrubland is a defensible first pass only if you say so.
  • Over-filling sinks in karst or genuinely closed-basin terrain, which erases real internal drainage.

QA and validation

Before the layer leaves your machine: confirm all five input rasters share CRS, extent, and cell size (gdalinfo on each); sanity-check that A is zero where C or LS is zero and never negative; histogram the result and confirm the high-loss tail aligns with visibly steep, bare slopes on the hillshade or imagery; and spot-check a handful of cells by hand as in the worked example. If you have any monitoring or field-survey points, compare ranked order rather than absolute tonnage, because RUSLE's absolute values carry large uncertainty.

Bathyl perspective

We treat an erosion surface as a prioritisation layer with a stated model, stated inputs, and stated limits, never as a stand-alone verdict. A map that says "RUSLE sheet-and-rill estimate, 10 m grid, C from 2024 WorldCover, validate gully reaches in the field" is far more useful to an engineer than a glossy red-orange raster with no provenance.

Related reading

Sources