What layers belong in a terrain analysis stack and in what order?

Start from a clean, projected DEM, then derive slope and aspect, hillshade for visualization, curvature and ruggedness (TRI/TPI) for landform character, and flow direction/accumulation for drainage. Build measurement layers from the DEM directly, not from each other.

Why must a DEM be in a projected CRS before computing slope?

Slope is rise over run, and run must be a real horizontal distance. If the DEM is in a geographic CRS, the horizontal unit is degrees while elevation is metres, so the ratio is meaningless. Reproject to a metric projected CRS first.

What is the z-factor and when do I need it?

The z-factor converts vertical units to the horizontal units used by the slope/hillshade math. It is 1 when both are metres, but must be set (for example about 0.3048 for feet elevation over metric ground, or a latitude-dependent value for geographic DEMs) when units differ.

How do I avoid noise and artifacts in terrain derivatives?

Match DEM resolution to the question rather than always picking the finest data, fill voids and remove spikes before deriving products, use an appropriate neighborhood, and inspect for striping, edge effects, and pits before trusting outputs.

Building a Terrain Layer Stack From a DEM

The short version

A terrain layer stack is the ordered set of products you derive from a digital elevation model — slope, aspect, hillshade, curvature, ruggedness, and drainage — built so each one is computed on a clean, correctly projected DEM with the right units. The single most important rule: derive measurement layers (slope, aspect, curvature) from the DEM directly, in a projected metric CRS, with a correct z-factor. Hillshade is the exception — it is for the eye, not the calculator, and should never be used as a measurement.

Most terrain mistakes are not in the algorithm; they are in the inputs. A DEM in the wrong CRS, with elevation in feet over ground in metres, with unfilled voids, or at a resolution that does not match the question, produces slope rasters that look plausible and are wrong. Get the DEM right first and the derivatives follow.

Start with the DEM, and check it

Before deriving anything, confirm five properties of the source DEM:

It is elevation data, not a rendered hillshade or coloured relief image. A pre-shaded GeoTIFF carries pixel brightness, not height; you cannot recover slope from it.
Horizontal CRS is projected and metric. Slope is rise/run; the run must be a real distance. A DEM in EPSG:4326 has a horizontal unit of degrees, so reproject:
```
gdalwarp -t_srs EPSG:25831 -r bilinear -tr 10 10 dem_4326.tif dem_25831.tif
```
Use bilinear (or cubic) resampling for continuous elevation, never nearest neighbour, which preserves blocky steps that become false terraces in slope.
Vertical units are known. If elevation is in feet and ground is in metres, you must set a z-factor (below).
NoData is declared and voids are handled. Undeclared NoData (e.g. -9999 treated as real) creates cliffs at the data edge.
Resolution matches the question. A 1 m lidar DEM resolves individual boulders and microtopography; for regional landform screening that detail is noise. Choosing the finest available DEM by reflex is a common error.

Clean before you derive

Two preprocessing steps prevent most artifacts:

Fill voids. Gaps in radar-derived DEMs (SRTM in steep terrain, for instance) produce spurious extremes. Use gdal_fillnodata or QGIS Fill nodata before deriving products.
Remove spikes/pits that are sensor noise, not real. For hydrological work, hydrological pit-filling (Wang & Liu, or the Fill sinks tools in SAGA/Whitebox) is required so flow does not pool in artifacts — but only apply it to the hydrology branch, because filling alters real closed depressions you might care about elsewhere.

The stack, layer by layer

Slope

Slope is the rate of elevation change. It can be expressed in degrees (0–90) or percent rise (rise/run × 100), and these are not interchangeable: 45° equals 100% rise, not 45%. State which you used.

gdaldem slope dem_25831.tif slope_deg.tif -compute_edges
gdaldem slope dem_25831.tif slope_pct.tif -p -compute_edges

-compute_edges avoids a NoData fringe at the raster border. In QGIS the equivalent is Raster Terrain Analysis > Slope; in ArcGIS Pro, Spatial Analyst Slope. All compute slope over a 3×3 neighbourhood (Horn's method by default in GDAL/ArcGIS), so cell size directly sets the scale of what you measure.

Aspect

Aspect is the compass direction the slope faces, 0–360° clockwise from north, with flat cells flagged (often -1 or -9999). It is circular data: never average aspect arithmetically (the mean of 350° and 10° is not 180°). Use it for solar exposure, hydrology, and habitat work.

gdaldem aspect dem_25831.tif aspect.tif -compute_edges

Hillshade (visualization only)

Hillshade simulates illumination from a light source. Defaults are azimuth 315° (NW) and altitude 45°. It is excellent for revealing landform structure and for basemaps under other layers, but it encodes lighting, not measurable terrain — do not treat a dark slope as a hazard value.

gdaldem hillshade dem_25831.tif hillshade.tif -az 315 -alt 45 -compute_edges

A multidirectional hillshade (-multidirectional) reduces the loss of detail on slopes facing away from a single light source.

The z-factor

Slope, aspect, and hillshade all assume vertical and horizontal units match. When they do not, set a z-factor that converts vertical to horizontal units. For feet elevation over metric ground, z ≈ 0.3048. For a geographic DEM left in degrees (not recommended for measurement), the horizontal degree length varies with latitude, so a single z-factor is only an approximation — another reason to reproject first. In gdaldem use -s (scale) for this; in ArcGIS/QGIS it is the explicit z-factor parameter.

Curvature and ruggedness

For landform character beyond slope:

Curvature (profile and plan) highlights ridges, valleys, and breaks of slope — useful for erosion and deposition reasoning.
TRI (Terrain Ruggedness Index) and TPI (Topographic Position Index) quantify local roughness and relative position. gdaldem TRI and gdaldem TPI compute them; TPI's meaning depends entirely on the neighbourhood size, so pick it to match the landform scale you care about.

Drainage

On the hydrologically-filled DEM, derive flow direction then flow accumulation (Whitebox, SAGA, or GRASS r.watershed). Accumulation thresholds give a synthetic stream network for screening. Keep this branch separate from the measurement branch because of the pit-filling difference noted above.

Worked example: a slope screen for a road corridor

Inputs: a 1 m lidar DEM in EPSG:4326, vertical in metres. The road question is about general grades over kilometres, not centimetre roughness.

Reproject and resample to a sensible scale: gdalwarp -t_srs EPSG:25831 -tr 5 5 -r cubic dem.tif dem_5m.tif. Going from 1 m to 5 m suppresses microtopographic noise that would otherwise dominate the slope histogram.
Fill voids: gdal_fillnodata dem_5m.tif dem_5m_filled.tif.
Slope in percent: gdaldem slope dem_5m_filled.tif slope_pct.tif -p -compute_edges.
Reclassify into grade bands (e.g. <8%, 8–15%, >15%) for screening.

Run on the raw 1 m DEM, the same screen flagged thousands of "steep" cells that were boulders and ditch edges — true detail, wrong scale for the question.

QA before you trust the stack

Compare hillshade against a known map of the area; landforms should read correctly.
Sanity-check slope against a hand calculation on a feature of known grade.
Inspect the raster border for NoData fringe (forgot -compute_edges) and the data edge for cliffs (undeclared NoData).
Look for striping or stepping in slope, which signals an integer DEM, wrong resampling, or sensor artifacts.
When comparing two derivatives, confirm identical cell size, extent, and method first.

Bathyl perspective

We build terrain stacks as evidence layers with explicit limits: each product records the source DEM, its resolution, CRS, units, and the method used, and we are clear about which layers measure (slope, aspect, curvature) and which only illustrate (hillshade). Matching DEM resolution to the actual question — rather than reaching for the finest data — is usually what separates a defensible terrain analysis from a noisy one.

Building a Terrain Layer Stack