Computing Slope and Aspect from a DEM with NumPy

Compute slope and aspect from a digital elevation model with NumPy and SciPy: the Horn kernel, correct cell-size scaling, nodata handling and circular aspect encoding.

Convolve the elevation array with the two Horn 3×3 kernels, divide each result by 8 * cell_size, then take degrees(arctan(hypot(dz_dx, dz_dy))) for slope and arctan2 for aspect. The cell size must be in metres, which means the DEM must be in a projected CRS before you start. Everything else — nodata, flat ground, edge handling — is bookkeeping, and getting the bookkeeping wrong is what produces slope maps that look plausible and are quietly off by a factor of a hundred.

This page covers the kernel mechanics in detail. For where slope and aspect sit in a wider feature stack — curvature, topographic position, roughness — see DEM and Terrain Derivative Features.

The two Horn kernels and what they measure Two three by three weight grids are shown side by side. The first has minus one, minus two, minus one in its left column and plus one, plus two, plus one in its right column, measuring the east-west gradient. The second has the same weights arranged by row, measuring the north-south gradient. Beneath each grid the corresponding division by eight times the cell size is written. Two convolutions, one divisor, and the units come out right kx — east–west gradient −1 0 +1 −2 0 +2 −1 0 +1 dz_dx = conv(dem, kx) / (8 × cell) weights sum to 8 per side ky — north–south gradient −1 −2 −1 0 0 0 +1 +2 +1 dz_dy = conv(dem, ky) / (8 × cell) slope = atan(hypot(dz_dx, dz_dy)) Drop the cell-size divisor and every slope is wrong by exactly that factor.

Why This Fails in Geospatial ML Pipelines

The first failure is units, and it is silent. Slope is a dimensionless ratio, so nothing in the code objects when the denominator is a fraction of a degree instead of a number of metres. A DEM in EPSG:4326 at roughly 0.0001° per pixel produces gradients scaled by about 1/11 of the correct value at the equator and by more nearer the poles — so the feature is not merely wrong, it is wrong in a way that varies systematically with latitude. A model fitted on that feature will learn latitude.

The second failure is nodata contamination. DEM tiles ship with sentinel values like -32768 or -9999. Convolution does not know they are special: a single sentinel cell drags the gradient of all eight of its neighbours towards a cliff, so the slope map grows a halo of implausible values around every void and along every tile seam. Because slope has no natural upper bound in the data, nothing flags it.

The third is the aspect wrap-around. Aspect is a bearing in [0, 360), and 359° is adjacent to 1° on the ground but 358 units away in the feature matrix. Tree models split on that gap; linear models read it as an enormous difference. The remedy — decomposing to northness and eastness — is two lines, and it typically matters more to model performance than the choice of estimator, which is the same reasoning applied to categorical variables in encoding categorical geographic features.

Core Principles

  • Validate the CRS before differentiating. Refuse to run on an angular CRS. The check is three lines and prevents the most expensive class of error here.
  • Convert nodata to NaN at read time. Sentinel values must leave the numeric domain before any kernel touches them, not be filtered afterwards.
  • Divide by 8 * cell_size, once. The eight normalises the Horn weights; the cell size converts to a true gradient. Both, in that order.
  • Treat flat ground as directionless. Mask aspect where slope is below a tolerance, then encode those cells as zero in both circular components.
  • Emit circular components, not degrees. Keep the bearing for maps and reports; give the model cos and sin.
  • Test against a synthetic plane. A plane with a known grade is the only test that catches a units error, and it runs in milliseconds.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
import rasterio
from scipy import ndimage

logger = logging.getLogger(__name__)

HORN_X = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype="float32")
HORN_Y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype="float32")
FLAT_TOL = 1e-6          # gradients below this are treated as flat


def read_dem_metric(path: str) -> tuple[np.ndarray, float]:
    """Read a DEM as float32 with nodata as NaN, and return its cell size in metres.

    Raises:
        ValueError: if the DEM has no CRS, is geographic, or has non-square pixels.
    """
    with rasterio.open(path) as src:
        if src.crs is None:
            raise ValueError(f"{path}: no CRS — cannot scale gradients")
        if not src.crs.is_projected:
            raise ValueError(
                f"{path}: CRS is geographic ({src.crs.to_string()}). Reproject to a "
                "metric CRS such as a UTM zone before computing slope."
            )
        dem = src.read(1, masked=True).filled(np.nan).astype("float32")
        xres, yres = abs(src.transform.a), abs(src.transform.e)

    if not np.isclose(xres, yres, rtol=1e-3):
        raise ValueError(f"{path}: non-square pixels {xres} x {yres}; resample first")
    logger.info("read %s: %s px at %.3f m", path, dem.shape, xres)
    return dem, float(xres)


def slope_aspect(dem: np.ndarray, cell_size: float) -> tuple[np.ndarray, np.ndarray]:
    """Horn slope (degrees) and aspect (degrees clockwise from north).

    Aspect is NaN wherever the surface is flat, because direction is undefined there.
    """
    if cell_size <= 0:
        raise ValueError(f"cell_size must be positive, got {cell_size}")

    dz_dx = ndimage.convolve(dem, HORN_X, mode="nearest") / (8.0 * cell_size)
    dz_dy = ndimage.convolve(dem, HORN_Y, mode="nearest") / (8.0 * cell_size)

    grade = np.hypot(dz_dx, dz_dy)
    slope = np.degrees(np.arctan(grade)).astype("float32")

    aspect = np.degrees(np.arctan2(dz_dy, -dz_dx))
    aspect = ((450.0 - aspect) % 360.0).astype("float32")
    aspect[grade < FLAT_TOL] = np.nan
    return slope, aspect


def circular_components(aspect_deg: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Model-safe aspect: (northness, eastness) in [-1, 1], zero where undefined."""
    rad = np.radians(aspect_deg)
    northness = np.cos(rad).astype("float32")
    eastness = np.sin(rad).astype("float32")
    undefined = np.isnan(aspect_deg)
    northness[undefined] = 0.0
    eastness[undefined] = 0.0
    return northness, eastness

Step-by-Step Walkthrough

Step 1 — Read the DEM and confirm it is metric

dem, cell = read_dem_metric("dem_utm33n.tif")
print(f"{np.isnan(dem).mean():.2%} nodata, cell size {cell} m")

If this raises, reproject first. Reprojection resamples elevation, so do it once and cache the result rather than repeating it per run — the reprojection mechanics are in reprojecting a raster with rasterio warp reproject.

Step 2 — Deal with voids before convolving

void = np.isnan(dem)
if void.any():
    idx = ndimage.distance_transform_edt(void, return_distances=False, return_indices=True)
    dem_filled = dem.copy()
    dem_filled[void] = dem[tuple(i[void] for i in idx)]
else:
    dem_filled = dem

Nearest-neighbour fill is deterministic and cannot invent elevations outside the observed range. Keep the original void mask — you will want to restore NaN in the outputs so the model sees missing rather than fabricated terrain.

Step 3 — Compute slope and aspect

slope, aspect = slope_aspect(dem_filled, cell)
slope[void] = np.nan
aspect[void] = np.nan
print(f"slope: {np.nanmin(slope):.1f}–{np.nanmax(slope):.1f}°, "
      f"median {np.nanmedian(slope):.1f}°")

A median slope of 0.02° over mountainous terrain, or a maximum of 89.99° over farmland, both point at the cell-size divisor.

Step 4 — Encode aspect for the model

northness, eastness = circular_components(aspect)
features = np.stack([slope, northness, eastness]).astype("float32")
assert features.shape[0] == 3

Step 5 — Write a three-band feature raster

with rasterio.open("dem_utm33n.tif") as src:
    profile = src.profile.copy()
profile.update(count=3, dtype="float32", nodata=np.nan,
               compress="deflate", tiled=True, blockxsize=512, blockysize=512)

with rasterio.open("slope_aspect.tif", "w", **profile) as dst:
    for i, (band, name) in enumerate(
            zip(features, ["slope_deg", "northness", "eastness"]), start=1):
        dst.write(band, i)
        dst.set_band_description(i, name)
One untreated nodata cell corrupts a whole neighbourhood A five by five grid of elevation values contains one cell holding the sentinel value minus nine thousand nine hundred and ninety-nine. The eight cells surrounding it are highlighted as contaminated, because the three by three kernel centred on each of them reads the sentinel and reports a near-vertical slope. A sentinel value is just a number to a convolution 412 414 417 419 421 411 413 416 418 420 410 412 −9999 417 419 409 411 414 416 418 8 contaminated neighbours each reports a near-vertical slope One void becomes a 3×3 patch of nonsense — and a tile seam becomes a wall of it. Read with masked=True and fill with NaN before any kernel runs.

Verification

A synthetic plane is the whole test suite for units, and it takes microseconds.

import numpy as np


def test_slope_matches_known_grade():
    """A 10% grade must give atan(0.1) = 5.711 degrees, whatever the cell size."""
    for cell in (1.0, 10.0, 30.0):
        cols = 40
        plane = (np.arange(cols, dtype="float32") * cell * 0.1)[None, :].repeat(40, 0)
        slope, aspect = slope_aspect(plane, cell)
        assert np.allclose(slope[2:-2, 2:-2], np.degrees(np.arctan(0.1)), atol=1e-3)
        # rising towards +x means the surface faces west: bearing 270
        assert np.allclose(aspect[2:-2, 2:-2], 270.0, atol=1e-3)


def test_flat_surface_has_undefined_aspect():
    flat = np.full((20, 20), 100.0, dtype="float32")
    slope, aspect = slope_aspect(flat, 10.0)
    assert np.allclose(slope, 0.0, atol=1e-6)
    assert np.isnan(aspect).all()
    n, e = circular_components(aspect)
    assert np.allclose(n, 0.0) and np.allclose(e, 0.0)


def test_matches_gdaldem(tmp_path):
    """Optional cross-check against the reference implementation."""
    import subprocess
    subprocess.run(["gdaldem", "slope", "dem_utm33n.tif", str(tmp_path / "ref.tif")],
                   check=True)
    with rasterio.open(tmp_path / "ref.tif") as src:
        ref = src.read(1, masked=True).filled(np.nan)
    ours, _ = slope_aspect(*read_dem_metric("dem_utm33n.tif"))
    inner = (slice(2, -2), slice(2, -2))
    assert np.nanmax(np.abs(ours[inner] - ref[inner])) < 0.01

If the plane test passes and real data still looks wrong, the problem is upstream: check for a vertical unit in feet, or a DEM that was resampled with an averaging filter that flattened the terrain.

The order these steps run in is not negotiable A left-to-right flow of six stages: validate that the CRS is projected, convert the nodata sentinel to NaN, fill small voids, convolve with the Horn kernels and divide by eight times the cell size, mask aspect on flat cells, and finally split aspect into northness and eastness for the model. Reordering any two of these produces a plausible, wrong feature 1. CRS projected? 2. nodata → NaN 3. fill voids keep the mask 4. convolve ÷ 8 × cell 5. mask flat cells 6. split cos / sin Fill before convolve, or voids spread. Mask after convolve, or flat cells claim north. Test with a synthetic plane after any change to steps 4 or 6.

FAQ

Why divide the convolution result by eight times the cell size?

The Horn weights sum to eight on each side of the centre, so dividing by eight recovers the mean elevation difference per neighbour; dividing by the cell size converts that difference into rise per unit run. Skip either divisor and slope is wrong by exactly that constant factor — which looks like a plausible map and destroys the feature.

What should aspect be on flat ground?

Undefined. arctan2(0, 0) returns zero, which would be reported as due north, so a plain would appear to face north everywhere. Mask aspect to NaN below a small gradient tolerance and encode those cells as zero northness and zero eastness.

Does NumPy give the same answer as gdaldem slope?

Yes, to floating point precision, as long as you use the Horn weights, divide by 8 * cell_size, and handle edges the same way. That makes gdaldem a convenient reference in tests, as in the third test above.

Part of: DEM and Terrain Derivative Features Part of: Spatial Feature Engineering for Machine Learning