Short answer

To compute slope from a digital elevation model with GDAL, run gdaldem slope dem.tif slope.tif. By default gdaldem returns slope in degrees (0–90), computed with Horn's 3×3 finite-difference algorithm over the eight neighbours of each cell. The two mistakes that wreck most slope rasters are running the command on a DEM in a geographic CRS (degrees on X/Y, metres on Z) and forgetting that the result inherits the DEM's vertical noise. Get the CRS and units right first, and the slope follows.

What slope actually measures

Slope is the magnitude of the maximum rate of elevation change at a cell. gdaldem estimates the partial derivatives dz/dx and dz/dy from a moving 3×3 window, then combines them. In degrees:

slope_deg = atan( sqrt( (dz/dx)^2 + (dz/dy)^2 ) ) * 180/pi

In percent it is simply 100 * sqrt((dz/dx)^2 + (dz/dy)^2). A 45° slope is 100% — a fact that surprises people who expect "100%" to mean vertical. Percent slope is unbounded and climbs fast: a 60° wall is about 173%, an 80° cliff about 567%. For engineering, hydrology, and land-capability work, degrees are usually safer because they are bounded and intuitive; percent is common in road design and erosion (USLE/RUSLE) contexts.

The default algorithm is Horn (1981), which weights the diagonal neighbours by 1 and the cardinal neighbours by 2. The alternative is Zevenberger–Thorne, selected with -alg ZevenbergenThorne, which fits a partial quartic surface and uses only the four cardinal neighbours. Horn is more resistant to noise (it smooths a little); Zevenberger–Thorne tracks fine structure more faithfully but amplifies DEM noise. For lidar-derived 1 m DEMs that are already crisp, Zevenberger–Thorne is defensible; for resampled or interpolated SRTM/Copernicus data, keep Horn.

The horizontal–vertical unit trap

This is the single most common slope error. gdaldem computes the gradient as dz divided by dx, where dx is the cell size taken from the geotransform. If your DEM is in EPSG:4326, the cell size is expressed in degrees of longitude/latitude, but dz is in metres. Dividing metres by degrees produces a number with no physical meaning, and your slope raster will be wildly wrong — typically far too shallow.

There are two correct fixes:

Fix 1 — reproject to a metric CRS (preferred). Pick a projected CRS appropriate to the area: a UTM zone (e.g. EPSG:32633 for UTM 33N), a national grid (EPSG:25832 ETRS89 / UTM 32N for central Europe, EPSG:2154 RGF93 / Lambert-93 for France), or an equal-area CRS for continental extents.

gdalwarp -t_srs EPSG:25832 -tr 30 30 -r bilinear \
  -of GTiff dem_4326.tif dem_utm.tif
gdaldem slope dem_utm.tif slope_deg.tif -alg Horn

Use bilinear (or cubic) resampling for continuous elevation, never near — nearest-neighbour creates stair-steps that show up as false slope spikes.

Fix 2 — use -scale. If you must keep a geographic CRS, tell gdaldem how many ground units correspond to one horizontal dataset unit:

gdaldem slope dem_4326.tif slope_deg.tif -scale 111120

111120 is the approximate number of metres in one degree of latitude (the WGS84 mean meridian length / 360). This is only exact in the north–south direction; longitude degrees shrink with cos(latitude), so -scale introduces a slope error that grows toward the poles. At 45° latitude, east–west distances are overstated by roughly 1/cos(45°) ≈ 1.41, biasing slope low on east–west facing terrain. For anything beyond a quick look, reproject. The GDAL docs spell this out under the gdaldem reference.

A complete worked workflow

Say you receive a 1 arc-second SRTM tile covering a mountainous study area and need a slope layer for a landslide-susceptibility screen.

  1. Inspect before touching it. Confirm CRS, resolution, NoData, and value range:
gdalinfo -stats srtm_n46e007.tif

Look at Coordinate System (expect EPSG:4326), Pixel Size, and NoData Value (SRTM voids are often -32768). If the NoData is not set in the header, fix it with gdal_edit.py -a_nodata -32768 srtm_n46e007.tif so voids do not contaminate the gradient.

  1. Fill voids if present. Unfilled voids create rings of extreme slope. Interpolate small holes:
gdal_fillnodata.py -md 10 srtm_n46e007.tif srtm_filled.tif
  1. Reproject to UTM at the native ground resolution (~30 m for 1 arc-second near the equator; choose the right UTM zone):
gdalwarp -t_srs EPSG:32632 -tr 30 30 -r cubic \
  -dstnodata -32768 srtm_filled.tif srtm_utm.tif
  1. Compute slope in degrees, and separately in percent if needed:
gdaldem slope srtm_utm.tif slope_deg.tif -compute_edges
gdaldem slope srtm_utm.tif slope_pct.tif -p -compute_edges

-compute_edges is important: without it the outermost row and column of pixels get NoData because the 3×3 window runs off the raster. With it, GDAL uses a reduced stencil at the border so you keep full coverage — valuable when tiles will be mosaicked.

  1. Reclassify into susceptibility bands if that is the goal. A common scheme buckets slope into stability classes (for example 0–5°, 5–15°, 15–25°, 25–35°, >35°). Use gdal_calc.py or a QGIS reclassify; the exact thresholds should come from your geotechnical model, not a default.

Common pitfalls and why they happen

  • Slope looks "too smooth" or "too noisy." The DEM resolution, not the slope command, controls this. Slope computed on a 30 m DEM averages microtopography; on a 1 m lidar DEM it captures boulders and ditches. Pick the DEM resolution that matches your phenomenon before blaming the algorithm.
  • Banding or terracing in the slope raster. This is quantisation: integer DEMs (e.g. metres stored as Int16) produce flat steps that turn into stripes of constant slope. Convert to Float32 (gdal_translate -ot Float32) before differencing, and resample with cubic not near.
  • A border of NoData around every tile. You forgot -compute_edges. Add it, or expect seams after mosaicking.
  • Slope of zero across an apparently varied area. Usually the unit trap (geographic CRS without -scale) collapsing real relief into near-flat values.
  • Extreme spikes along rivers or coastlines. Often unfilled NoData voids or a hydro-flattening artifact in the source DEM, not real terrain.

Validation

Sanity-check the output statistics:

gdalinfo -stats slope_deg.tif

For a degree slope, Minimum should be 0 and Maximum must not exceed 90 — a max above 90 is a guaranteed sign of a unit or NoData problem. Spot-check a few cells against a manual estimate: pick two points a known distance apart on a uniform hillside, read their elevations, and compute atan(rise/run) by hand; it should land within a degree or two of the raster. Finally, overlay the slope raster on a hillshade — slope highs should hug the steep flanks visible in the shading, and flat valley floors and ridge crests should read near zero.

Bathyl perspective

A slope layer is only as trustworthy as the DEM and the CRS beneath it. We treat the reprojection, void handling, and unit check as part of the deliverable, not preamble, and we record the exact gdaldem invocation alongside the output so any downstream susceptibility or routing model can be reproduced from the raw tile.

Related reading

Sources