Generate Slope, Aspect and Hillshade in PyQGIS

Slope, aspect and hillshade come from the same three-by-three neighbourhood calculation, and all three are a single algorithm call. What separates a usable result from a convincing-looking wrong one is entirely in the setup: the CRS the DEM sits in, the vertical unit it stores, and — for hillshade only — a couple of cartographic conventions that readers' eyes depend on.

This recipe belongs to Terrain & Interpolation Analysis in PyQGIS. It covers producing each raster from Python, choosing between the native and GDAL implementations, handling flat ground in an aspect raster, and building a multi-directional hillshade that shows features a single light source hides.

Three derivatives, one set of prerequisitesBefore deriving anything the DEM must be in a projected CRS with vertical units matching, and its nodata must be declared. From there slope produces degrees or percent, aspect produces a bearing that is undefined on flat ground, and hillshade produces an eight bit image intended for display rather than measurement.Two of these are measurements; one is a picturecheck firstprojected CRS?vertical unit = horizontal?nodata declared?skip these and everything below is wrongslopeFloat32 · degrees or percentaspectFloat32 · 0–360, circularhillshadeByte · 0–255, for displaythresholds, suitability,zonal summariessolar exposure, habitat,reclassified to compass pointsunderneath a thematic fillset to multiply

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer, with the GDAL provider available (it is, in every standard install).
  • A DEM in a projected CRS whose horizontal units match its vertical units. Reproject first if not — see batch reprojecting raster datasets.
  • Nodata declared on the input. An undeclared fill value of -9999 becomes a cliff face several kilometres high at every gap.

Validate before you derive

Two lines of checking save an afternoon of confused output.

from qgis.core import QgsRasterLayer

dem = QgsRasterLayer("/data/dem_27700.tif", "dem")
if not dem.isValid():
    raise RuntimeError("DEM failed to load — check the path and driver")
if dem.crs().isGeographic():
    raise RuntimeError(f"{dem.crs().authid()} is geographic; reproject before deriving")
if not dem.dataProvider().sourceHasNoDataValue(1):
    dem.dataProvider().setNoDataValue(1, -9999)

Breakdown: isValid() catches the case where the path is wrong or the driver is missing, which otherwise surfaces much later as an empty output rather than an error. The geographic check is the one that matters most: in EPSG:4326 the cell size is in degrees and the heights are in metres, so slope comes out around five orders of magnitude too small — a hillside reported as 0.0003 degrees. Declaring nodata is done on the layer rather than the file, which is enough for the Processing algorithms that receive the layer object.

Passing the layer object rather than a path is generally the better habit here, because the nodata declaration and any other layer-level settings travel with it.

Slope, in degrees or percent

The native algorithm returns degrees; GDAL's returns degrees or percent.

import processing

processing.run("native:slope", {
    "INPUT": dem,
    "Z_FACTOR": 1.0,
    "OUTPUT": "/data/output/slope_deg.tif",
})

processing.run("gdal:slope", {
    "INPUT": dem,
    "BAND": 1,
    "SCALE": 1.0,
    "AS_PERCENT": True,
    "COMPUTE_EDGES": False,
    "ZEVENBERGEN": False,
    "OUTPUT": "/data/output/slope_pct.tif",
})

Breakdown: GDAL calls the vertical scaling SCALE where the native algorithm calls it Z_FACTOR; both mean "multiply the heights by this before differentiating", and both should be 1.0 unless the DEM's vertical unit differs from the CRS's horizontal one — 0.3048 for a DEM in feet on a metric grid. ZEVENBERGEN switches from the default Horn algorithm to Zevenbergen-Thorne, which weights the diagonal neighbours differently and gives slightly sharper results on smooth surfaces; Horn is more robust on noisy DEMs and is the right default. AS_PERCENT matters because engineering criteria are usually written as percentages: a 1-in-4 gradient is 25 percent and 14 degrees, and confusing the two is a real and expensive mistake.

Aspect, and the flat-ground problem

Aspect is the compass bearing of the downhill direction, and on perfectly flat ground there is no such direction.

processing.run("native:aspect", {
    "INPUT": dem,
    "Z_FACTOR": 1.0,
    "OUTPUT": "/data/output/aspect.tif",
})

Breakdown: Flat cells receive -1 in GDAL's convention and -9999 in some others, so any reclassification must handle the flag explicitly rather than letting it fall into the "north" bin — where it will, since a naive classification of 0–360 into eight sectors puts everything below 22.5 degrees into north and a -1 sorts below that.

The second trap is that aspect is circular. A mean of 350 and 10 degrees is 0, not 180, so averaging an aspect raster with zonal statistics produces meaningless numbers. When a summary is needed, reclassify into compass sectors first and take the majority:

processing.run("native:reclassifybytable", {
    "INPUT_RASTER": "/data/output/aspect.tif",
    "RASTER_BAND": 1,
    "TABLE": [
        0, 22.5, 1,  22.5, 67.5, 2,  67.5, 112.5, 3,  112.5, 157.5, 4,
        157.5, 202.5, 5,  202.5, 247.5, 6,  247.5, 292.5, 7,
        292.5, 337.5, 8,  337.5, 360, 1,
    ],
    "NO_DATA": -9999,
    "RANGE_BOUNDARIES": 0,
    "OUTPUT": "/data/output/aspect_class.tif",
})

Breakdown: The table is a flat list of min, max, value triples, and north appears twice because it wraps around zero — that wrap is the whole reason this reclassification exists. RANGE_BOUNDARIES: 0 makes each range include its minimum and exclude its maximum, which avoids a cell exactly on a boundary landing in two classes. With classes in hand, the majority statistic in zonal statistics gives a defensible "this catchment mostly faces south-west".

Aspect wraps, and flat ground is not a directionThe eight compass sectors each span forty five degrees. North is split across the wrap at zero, so it must appear twice in any reclassification table. Cells on flat ground carry a sentinel value rather than a bearing and belong in their own class, not in the north sector.North is one class written as two rangesNESW337.5–360plus 0–22.5eight sectors, nine rangesthe extra range is north, past the wrapboth map to class 1flat cellsvalue −1, not a bearingexclude, or they all become north

Hillshade, and why azimuth 315 is not negotiable

Hillshade is a rendering of a light source over the surface. Two of its parameters are cartographic conventions rather than choices.

processing.run("native:hillshade", {
    "INPUT": dem,
    "Z_FACTOR": 2.0,
    "AZIMUTH": 315,
    "V_ANGLE": 45,
    "OUTPUT": "/data/output/hillshade.tif",
})

Breakdown: An azimuth of 315 degrees puts the light in the north-west. Human perception assumes light comes from above and to the left; lighting from the south-east makes most readers see ridges as valleys and valleys as ridges, an illusion strong enough that it does not go away once you know about it. V_ANGLE of 45 degrees is a balanced default — lower angles lengthen shadows and reveal subtle features at the cost of saturating steep ground to black. A Z_FACTOR above 1 exaggerates relief, which is legitimate here precisely because a hillshade is not a measurement.

Where the terrain has strong linear structure — glacial valleys, dune fields, fault scarps — a single light direction hides everything parallel to it. Combining several directions fixes that.

import processing

shades = []
for azimuth in (315, 15, 75, 135):
    shades.append(processing.run("native:hillshade", {
        "INPUT": dem, "Z_FACTOR": 2.0, "AZIMUTH": azimuth, "V_ANGLE": 45,
        "OUTPUT": "TEMPORARY_OUTPUT",
    })["OUTPUT"])

processing.run("native:cellstatistics", {
    "INPUT": shades,
    "STATISTIC": 2,                    # mean
    "OUTPUT": "/data/output/hillshade_multi.tif",
})

Breakdown: native:cellstatistics combines aligned rasters cell by cell, and the mean of four hillshades is the standard multi-directional recipe — softer than any single one, with no direction of structure invisible. The inputs are temporary outputs, so nothing but the final file is written to disk. Because all four derive from the same DEM they are automatically aligned, which is the condition cellstatistics requires and does not check.

Styling the outputs so they read correctly

Each of the three derivatives wants a different treatment, and the defaults QGIS applies are wrong for two of them.

Slope suits a sequential ramp with class breaks at the thresholds that matter to the work rather than at equal intervals — 0–2, 2–5, 5–10, 10–20, above 20 degrees says something about buildability, while five equal classes says only that the algorithm ran. Aspect needs a cyclic ramp, because a sequential one puts a hard colour break at north where the values wrap and invents a boundary that is not in the terrain. Hillshade wants greyscale with no ramp at all, sitting underneath everything else.

from qgis.core import (
    QgsRasterLayer, QgsProject, QgsSingleBandGrayRenderer, QgsContrastEnhancement,
)
from qgis.PyQt.QtGui import QPainter

shade = QgsRasterLayer("/data/output/hillshade_multi.tif", "hillshade")
QgsProject.instance().addMapLayer(shade)

renderer = QgsSingleBandGrayRenderer(shade.dataProvider(), 1)
enhancement = QgsContrastEnhancement(shade.dataProvider().dataType(1))
enhancement.setContrastEnhancementAlgorithm(
    QgsContrastEnhancement.StretchToMinimumMaximum
)
enhancement.setMinimumValue(0)
enhancement.setMaximumValue(255)
renderer.setContrastEnhancement(enhancement)
shade.setRenderer(renderer)
shade.triggerRepaint()

Breakdown: Pinning the stretch to the full 0–255 range rather than letting QGIS compute it from the visible extent is what stops the hillshade re-contrasting itself as you pan, which otherwise makes flat regions look dramatic and mountains look washed out. The layer belongs at the bottom of the tree with the thematic layers above it set to multiply, as covered in opacity and blend modes — that arrangement, rather than any transparency setting, is what produces shaded thematic mapping.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.16 LTR3.7native:slope, native:aspect, native:hillshade and the GDAL equivalents present.
3.22 LTR3.9native:cellstatistics available for combining aligned rasters.
3.28 LTR3.9native:reclassifybytable boundary handling stable.
3.34 LTR3.12Baseline for this page.
3.40+3.12Additional terrain algorithms exposed through the native provider.

Troubleshooting

  • Slope is almost all zero. The DEM is in a geographic CRS. Reproject with bilinear resampling.
  • Slope values are ten times too large. Vertical and horizontal units differ. Set Z_FACTOR to the conversion factor, not to taste.
  • A cliff appears around every gap. Nodata was never declared, so the fill value is being treated as elevation.
  • The relief looks inverted. The azimuth is in the southern half. Return it to 315.
  • Averaging aspect gives nonsense. Aspect is circular. Reclassify into sectors and take the majority instead.
  • Everything flat is classified as north. The -1 flat sentinel fell into the first range. Exclude it explicitly in the reclassification table.
  • cellstatistics refuses the inputs. The rasters are not aligned. Derive them all from one DEM, or warp to a common grid first.

Conclusion

Check the CRS and the nodata before anything else, keep Z_FACTOR as a unit conversion rather than a style control, treat aspect as circular data with a separate flat class, and reserve exaggeration and multi-directional tricks for the hillshade, which is the only one of the three that is a picture rather than a measurement.

Frequently Asked Questions

Which is better, the native or GDAL slope? They agree to within rounding. GDAL is faster on large files and offers percent output, edge computation and the Zevenbergen-Thorne algorithm; native is simpler and returns a layer that slots into a chained workflow with less ceremony.

How do I get slope in a ratio like 1 in 20? Compute percent slope and take the reciprocal in the raster calculator: 100 / "slope_pct@1" gives the run for one unit of rise, with a guard needed where slope is zero.

Should I smooth the DEM before deriving slope? Only if the DEM is noisy and the noise is not real terrain. A light low-pass filter removes stepping artefacts from integer-valued DEMs; smoothing a good DEM discards real detail and flatters the result.

Can I compute these on a virtual raster? Yes. A VRT mosaic behaves as an ordinary input, which is the usual way to derive across a directory of tiles without building a physical mosaic first.