Python LiDAR

Generating Slope and Aspect Rasters with gdaldem

TL;DR: gdaldem slope dtm.tif slope.tif -alg Horn -compute_edges and gdaldem aspect dtm.tif aspect.tif -zero_for_flat, both on a DTM whose vertical and horizontal units match — and remember that slope is a magnitude you can average and aspect is a direction you cannot.

# Context and Motivation

This guide is part of Hillshade, Slope and Aspect, which covers the derivatives a terrain model supports. Hillshade is for looking at; slope and aspect are for computing with, and that difference is where the mistakes live.

Both come from the same two gradients that hillshade uses — the Horn stencil over a three by three neighbourhood. Slope is the magnitude of that gradient, in degrees or percent. Aspect is its direction, in degrees clockwise from north. They are produced by one pass each over the raster and cost almost nothing. What costs is using them without noticing that one of them is circular.

One gradient, two very different numbers The same terrain gradient produces two outputs. Slope is its magnitude, a value from zero to ninety degrees that can be averaged, differenced and compared like any other measurement. Aspect is its direction, a value from zero to 360 that wraps, so an ordinary mean of 350 and 10 degrees gives 180 — the exact opposite of the right answer. the Horn gradient dz/dx and dz/dy slope — a magnitude 0–90°, averages normally aspect — a direction 0–360°, wraps at north mean, deviation, difference between epochs — all valid mean(350°, 10°) = 180°, south when it should be north for aspect, average the unit vectors — atan2(mean(sin θ), mean(cos θ)) — or reclassify to sectors first and never feed a raw aspect raster to a tool that will resample it.

# Prerequisites and Assumptions

Requirement Detail
GDAL 3.x, for gdaldem
A DTM bare earth; slope over a DSM measures roof pitch
Matching units vertical and horizontal in the same unit, or a -s scale factor
A void-aware workflow NoData cells propagate to their neighbours in the derivatives

# Step-by-Step Implementation

# Step 1 — Confirm the units agree

A DTM in US survey feet on a metric grid needs -s 0.3048. Without it every slope is 3.28 times too steep — the failure shown in exporting hillshade from a LiDAR DTM.

# Step 2 — Compute slope

bash
gdaldem slope dtm.tif slope_deg.tif -alg Horn -compute_edges \
  -co COMPRESS=DEFLATE -co TILED=YES

Add -p for percent rise instead of degrees. -compute_edges fills the one-pixel border that would otherwise be NoData.

# Step 3 — Compute aspect

bash
gdaldem aspect dtm.tif aspect_deg.tif -compute_edges -zero_for_flat \
  -co COMPRESS=DEFLATE -co TILED=YES

Without -zero_for_flat, flat cells get −9999, which is correct and awkward; with it they get 0, which is north and wrong. Choose deliberately and record which you chose.

# Step 4 — Reclassify aspect before any statistics

Eight sectors, or a north-south exposure index, both behave arithmetically. Raw degrees do not.

# Complete Working Example

python
"""Derive slope and aspect, then summarise them correctly."""
from __future__ import annotations

import json
import subprocess
from pathlib import Path

import numpy as np
import rasterio


def derive(dtm: Path, out_dir: Path, z_factor: float = 1.0) -> dict[str, Path]:
    out_dir.mkdir(parents=True, exist_ok=True)
    slope = out_dir / "slope_deg.tif"
    aspect = out_dir / "aspect_deg.tif"
    common = ["-compute_edges", "-co", "COMPRESS=DEFLATE", "-co", "TILED=YES"]

    subprocess.run(["gdaldem", "slope", str(dtm), str(slope),
                    "-alg", "Horn", "-s", str(1.0 / z_factor), *common], check=True)
    subprocess.run(["gdaldem", "aspect", str(dtm), str(aspect),
                    "-zero_for_flat", *common], check=True)
    return {"slope": slope, "aspect": aspect}


def summarise(slope: Path, aspect: Path) -> dict:
    with rasterio.open(slope) as s:
        sl = s.read(1).astype("float64")
        sl_nodata = s.nodata
    with rasterio.open(aspect) as a:
        asp = a.read(1).astype("float64")

    valid = np.isfinite(sl)
    if sl_nodata is not None:
        valid &= sl != sl_nodata

    # Slope is a magnitude: ordinary statistics are fine.
    slope_stats = {
        "mean_deg": round(float(sl[valid].mean()), 2),
        "p95_deg": round(float(np.percentile(sl[valid], 95)), 2),
    }

    # Aspect is circular: average the unit vectors, not the degrees.
    sloped = valid & (sl > 2.0)  # a flat cell has no meaningful aspect
    theta = np.radians(asp[sloped])
    mean_dir = np.degrees(np.arctan2(np.sin(theta).mean(), np.cos(theta).mean())) % 360.0
    resultant = float(np.hypot(np.sin(theta).mean(), np.cos(theta).mean()))

    return {
        "slope": slope_stats,
        "aspect": {
            "mean_direction_deg": round(float(mean_dir), 1),
            "concentration": round(resultant, 3),
            "naive_mean_deg": round(float(asp[sloped].mean()), 1),
        },
    }


if __name__ == "__main__":
    paths = derive(Path("dtm.tif"), Path("derivatives"))
    print(json.dumps(summarise(paths["slope"], paths["aspect"]), indent=2))

The naive_mean_deg field is deliberate: printing both makes the difference visible to whoever reads the output next.

# Key Parameter Table

Flag Applies to Effect
-alg Horn slope, aspect The standard 3×3 estimator; ZevenbergenThorne is smoother
-p slope Percent rise instead of degrees
-s slope, hillshade Ratio of vertical to horizontal units; 0.3048 for feet over metres
-compute_edges all Fills the one-pixel NoData border
-zero_for_flat aspect Flat cells become 0 rather than −9999
-co all GeoTIFF creation options; compress and tile by default
Degrees and percent are not interchangeable Four gradients in both units. A gentle 5 degree slope is 8.7 percent. A moderate 15 degrees is 26.8 percent. A steep 30 degrees is 57.7 percent. And 45 degrees is exactly 100 percent, above which percent rise climbs without bound while degrees cannot exceed 90 — which is why a threshold written in one unit cannot be reused in the other. 5° — gentle 8.7% rise 15° — moderate 26.8% rise 30° — steep 57.7% rise 45° exactly 100% rise percent rise is unbounded above 45°, degrees stop at 90 — a ported threshold is always wrong

# Verification

Slope has no impossible values. Nothing above 90 degrees, and a maximum near 90 over natural terrain means the DTM has a spike.

Flat areas read as flat. Sample a car park or a lake margin; slope should be under a degree.

The circular mean differs from the naive one. On terrain with a genuine dominant aspect the two will disagree, and that disagreement is the whole reason for the extra code.

The derivative is of the DTM, not the DSM. Slope over a surface model measures roof pitch and canopy, which is occasionally what you want and usually not.

Why the naive mean lands in the wrong place A compass with a cluster of aspect values gathered around north, spanning 340 to 20 degrees. The circular mean sits at 358 degrees, inside the cluster. The naive arithmetic mean sits at 176 degrees, almost due south, because averaging numbers that wrap at 360 pulls the result to the middle of the numeric range rather than the middle of the directions. N S E W circular mean 358° — inside the cluster naive mean 176° — almost due south the same five values: 340°, 350°, 358°, 6°, 18° every one faces north; their arithmetic mean faces away from all of them.

# Gotchas and Edge Cases

Resampling an aspect raster produces nonsense. Interpolating between 355 and 5 gives 180. Resample the DTM and re-derive instead.

NoData spreads. A single void becomes a 3×3 block of NoData in every derivative. Fill voids in the DTM first — see filling NoData voids.

Slope depends on cell size. A 0.5 m DTM yields steeper slopes than a 2 m DTM of the same hill. Comparisons must use one grid.

Flat cells have no aspect. Mask by slope before using aspect for anything.

# Frequently Asked Questions

Why can I not average aspect values?

Because aspect wraps at 360 degrees. Five cells facing 340, 350, 358, 6 and 18 degrees all face north, and their arithmetic mean is 176 degrees — almost due south. Average the unit vectors instead, with atan2 of the mean sine over the mean cosine, or reclassify to sectors first.

Should I use -zero_for_flat?

It depends what consumes the raster, and the important thing is to decide rather than inherit a default. Without it flat cells are NoData, which is honest and awkward. With it they are zero, which is a legitimate aspect value meaning north, so a downstream statistic will treat every flat cell as north-facing.

Why does my slope raster look far too steep?

Almost always a unit mismatch: elevations in US survey feet on a metric grid, with no scale factor. Every slope then comes out 3.28 times too steep. Pass the ratio with -s, or better, rasterize in a CRS whose vertical unit matches its horizontal one.

Can I resample an aspect raster to a coarser grid?

No. Interpolating between 355 and 5 degrees gives 180, so resampling manufactures south-facing cells out of north-facing ones. Resample the DTM and re-derive the aspect from it.

Related in Hillshade, Slope and Aspect from LiDAR DTMs

↑ Hillshade, Slope and Aspect from LiDAR DTMs overview --- title: "Hillshade, Slope and Aspect from LiDAR DTMs" description: "Deriving hillshade, slope, and aspect rasters from a LiDAR-derived... Exporting Hillshade from a LiDAR DTM --- title: "Exporting Hillshade from a LiDAR DTM" description: "Step-by-step guide to producing a shaded-relief hillshade GeoTIFF from a...