Computing Slope and Aspect Rasters in Python

A slope raster returning values in the thousands of degrees, or a stability analysis flagging a flat field as a hazard, comes from one of three things: horizontal and vertical units that do not match, a window too small for the DEM’s noise, or NoData treated as an elevation. All three produce a raster that opens, renders and is wrong.

This page computes slope and aspect correctly, and covers the choices that decide whether the output is a measurement or an artefact.

Slope is a ratio, so the units must agree

Slope at a cell is the arctangent of the gradient magnitude: the rate of elevation change per unit horizontal distance. That ratio is only meaningful if the numerator and the denominator are in the same units.

In a projected CRS both are metres and the ratio is dimensionless. In a geographic CRS the horizontal spacing is in degrees while elevation is in metres, so the computed gradient is metres per degree — about 111,000 times too large near the equator, and varying with latitude. The result is not slightly wrong; it is wrong by five orders of magnitude, and by a different factor at each latitude.

The correct response is to reproject the DEM before differentiating. A vertical-exaggeration factor is sometimes offered as a workaround, and it is only correct at one latitude, so it converts an obvious error into a subtle one.

Why a geographic CRS breaks a slope computation Two cells side by side across a slope, shown in two coordinate reference systems. In a projected system both the horizontal spacing and the elevation difference are in metres, so the gradient is a dimensionless ratio and the slope is about eleven degrees. In a geographic system the horizontal spacing is expressed in degrees of longitude while the elevation difference is still in metres, so the computed ratio is enormous and the reported slope saturates near ninety degrees. A note adds that the conversion factor between the two depends on latitude, so no single vertical exaggeration corrects it across a survey. projected CRS — metres both ways Δx = 5.0 m Δz = 1.0 m slope = atan(1.0 / 5.0) = 11.3° geographic CRS — degrees and metres Δx = 0.000045° Δz = 1.0 m slope = atan(22 000) ≈ 90° The conversion factor varies with latitude, so a single vertical exaggeration is correct at one parallel and wrong elsewhere.

Figure 1 — The units failure. It is arithmetically obvious once seen and produces a raster that looks like a plausible near-vertical landscape, which is why it survives a quick glance.

Minimal reproducible solution

Refuse a geographic CRS outright, then differentiate over a window chosen from the DEM’s noise rather than from habit.

import numpy as np
import rasterio


def slope_aspect(dem_path: str, slope_path: str, aspect_path: str,
                 window_cells: int = 1) -> None:
    """Slope in degrees and aspect in degrees clockwise from north.

    window_cells > 1 differentiates over a wider stencil, which suppresses
    per-cell noise at the cost of resolving smaller features.
    """
    with rasterio.open(dem_path) as ds:
        if not ds.crs or not ds.crs.is_projected:
            raise ValueError(
                f"DEM is in {ds.crs}; slope needs matching horizontal and "
                "vertical units — reproject to a projected CRS first")
        z = ds.read(1, masked=True).filled(np.nan).astype(np.float64)
        px, py = abs(ds.transform.a), abs(ds.transform.e)
        profile = ds.profile

    k = window_cells
    # Central differences over a k-cell stencil; NaN propagates so NoData
    # neighbourhoods produce NoData rather than a fabricated gradient.
    dzdx = (np.roll(z, -k, axis=1) - np.roll(z, k, axis=1)) / (2 * k * px)
    dzdy = (np.roll(z, k, axis=0) - np.roll(z, -k, axis=0)) / (2 * k * py)
    dzdx[:, :k] = dzdx[:, -k:] = np.nan          # edges have no stencil
    dzdy[:k, :] = dzdy[-k:, :] = np.nan

    slope_deg = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
    aspect_deg = (np.degrees(np.arctan2(dzdy, -dzdx)) + 360.0) % 360.0
    aspect_deg[np.hypot(dzdx, dzdy) < 1e-9] = -1.0      # flat: aspect undefined

    profile.update(dtype="float32", nodata=np.nan, compress="deflate", predictor=3)
    for arr, path in ((slope_deg, slope_path), (aspect_deg, aspect_path)):
        with rasterio.open(path, "w", **profile) as dst:
            dst.write(arr.astype("float32"), 1)
            dst.update_tags(WINDOW_CELLS=str(window_cells))

Two details are easy to get wrong and hard to notice afterwards. Aspect is undefined on a flat cell — a zero gradient has no direction — and returning zero there means “north”, so every flat area reports as north-facing and any aspect-based analysis is biased. Marking it with a sentinel is the honest handling.

And NaN must propagate. Filling NoData with a number before differentiating manufactures an enormous gradient at every void edge, producing a ring of near-vertical slope around each hole. Letting NaN spread by one cell is correct: the gradient there genuinely is unknown.

Edge-case matrix

Input variant Naive result Correct handling
Geographic CRS Slope near 90° everywhere Refuse; reproject first
NoData filled with zero Ring of vertical slope at voids Propagate NaN
Flat terrain Aspect returns 0° = north Sentinel for undefined aspect
Noisy DEM, 3×3 window Slope noise of tens of degrees Widen the window, or smooth
Anisotropic cells Wrong gradient direction Use separate x and y spacings
Elevations in feet Slope wrong by 3.28× Convert before differentiating
Raster edges Wrap-around from np.roll Mask the border explicitly
Very steep terrain Fine Slope saturates gracefully at 90°

The anisotropic-cell row is worth a check rather than an assumption: a DEM whose pixels are not square — common after certain reprojections — needs its two spacings used separately, and code that uses a single cellsize returns an aspect rotated by an amount that varies with slope direction.

Verify the fix worked

Slope has a property that makes verification easy: it is bounded, and its distribution over a real landscape is characteristic.

import numpy as np
import rasterio


def assert_slope_plausible(slope_path: str, max_expected_deg: float = 75.0) -> None:
    """Bounds and distribution checks that catch every unit failure."""
    with rasterio.open(slope_path) as ds:
        s = ds.read(1, masked=True).compressed()

    assert s.size, "no valid slope values"
    assert float(s.min()) >= 0.0, "negative slope is impossible"
    assert float(s.max()) <= 90.0, "slope above 90° — check the arctangent"

    p99 = float(np.percentile(s, 99))
    assert p99 <= max_expected_deg, (
        f"99th percentile slope is {p99:.0f}° — a real landscape rarely has "
        "1% of its area steeper than this; suspect a unit mismatch")

    near_vertical = float((s > 85.0).mean())
    assert near_vertical < 0.02, (
        f"{near_vertical:.1%} of cells are near-vertical — either the DEM is "
        "in a geographic CRS or voids were filled before differentiating")


def assert_aspect_unbiased(aspect_path: str, tol: float = 0.35) -> None:
    """Aspect should not pile up on one compass direction."""
    with rasterio.open(aspect_path) as ds:
        a = ds.read(1, masked=True).compressed()
    a = a[a >= 0]                     # drop the undefined-aspect sentinel
    hist, _ = np.histogram(a, bins=8, range=(0, 360))
    frac = hist / max(hist.sum(), 1)
    assert frac.max() < 0.125 * (1 + tol) + 0.05, (
        f"aspect concentrates {frac.max():.1%} in one octant — flat cells are "
        "probably being reported as north-facing")

The aspect-bias test is the one that catches the flat-cell mistake, and nothing else does: a raster where every flat area reads as north has a perfectly valid range, a plausible mean, and a histogram with one enormous spike.

Aspect distribution with and without the flat-cell sentinel Two circular histograms of aspect over the same survey. In the first, flat cells were assigned zero degrees, so the north octant carries nearly half of all cells while the other seven share the remainder, producing an obviously biased rose. In the second, flat cells carry an undefined sentinel and are excluded, so the eight octants are close to even with a mild preference reflecting the site's real terrain. A note observes that only the shape of the distribution reveals the fault, because both rasters have identical valid ranges and plausible means. flat cells reported as north 47% of cells in one octant flat cells excluded roughly even, with real terrain preference Both rasters have a valid 0–360 range and a plausible mean; only the shape reveals the fault.

Figure 2 — The bias the range check cannot see. Aspect is circular, so a spike at one bearing is invisible in every summary statistic except the histogram.

When to escalate

  • Slope is plausible and a stability analysis disagrees with the ground. Slope from a DSM includes vegetation and structures; a stability analysis needs the bare-earth DTM. Confirm which surface was differentiated before questioning the numbers.
  • Widening the window enough to suppress noise erases the features being studied. The DEM’s noise is comparable to the relief of interest. Improving the DEM — denser imagery, a finer point cloud — is the only real fix; a wider window trades one problem for another.
  • Aspect is noisy on gentle slopes and clean on steep ones. That is inherent: aspect is the direction of a small vector, so its uncertainty grows as the gradient shrinks. Weight or mask aspect by slope magnitude rather than trying to stabilise it.

A note on what to deliver. Slope is far more useful to a client as a classified raster than as a continuous one, because the decisions it feeds are threshold-based: machinery access below a certain gradient, revegetation above another, batter angles in a specific band. Supplying both — the continuous float raster for anyone who wants to recompute, and a classified version whose breaks match the client’s own thresholds — costs one extra file and removes an entire round of “can you reclassify this for us” requests. Record the class breaks in the raster’s metadata so the classification is self-describing rather than depending on an email.

← Deriving Contours and Hillshade with GDAL

Slope noise and feature resolution against the differencing window Two curves against the half-width of the differencing stencil in cells. Slope noise, expressed as the standard deviation of computed slope over genuinely flat ground, falls steeply from about thirty degrees at a one-cell stencil to a few degrees by three cells. The smallest resolvable terrain feature grows linearly with the stencil, from about two cells to eight. A shaded band around two cells marks the region where noise is acceptable and features of interest still survive, with the note that the correct value depends on the DEM's own noise rather than on convention. workable 1 2 3 4 stencil half-width, in cells slope noise on flat ground smallest resolvable feature Measure the first curve on genuinely flat ground in your own DEM; the correct stencil follows from it, not from convention.

Figure 3 — The window is a trade between noise and resolution, and the position of the curves depends on the DEM. Measuring slope noise over a flat area is a one-line calibration that settles it.