Computing Slope and Aspect from a DEM

Use Horn’s weighted 3×3 differences, then convert gradient to slope and direction:

import numpy as np
from scipy import ndimage

KX = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) / 8.0   # east gradient
KY = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]]) / 8.0   # north gradient (row 0 is north)

dzdx = ndimage.convolve(dem, KX, mode="nearest") / res
dzdy = ndimage.convolve(dem, KY, mode="nearest") / res
slope = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
aspect = (np.degrees(np.arctan2(-dzdx, -dzdy)) + 360) % 360   # downslope, clockwise from north

The kernel signs and the aspect formula are where most implementations go wrong, and a wrong sign produces a plausible-looking but inverted result. This page belongs to terrain analysis and DEM-derived products in Satellite Processing Workflows & Index Pipelines.

The 3×3 Neighbourhood

Horn's weights For the east gradient, the left column of the neighbourhood is weighted minus one, minus two, minus one and the right column plus one, plus two, plus one, with the centre column ignored. For the north gradient the same weights run down the rows, with the top row positive because row zero is northernmost. The centre-adjacent cells count double, which suppresses noise while keeping sensitivity. east gradient weights north gradient weights -1 0 +1 -2 0 +2 -1 0 +1 +1 +2 +1 0 0 0 -1 -2 -1 Divide each sum by 8 × cell size. Top row positive because row 0 is northernmost in a north-up raster. Flip that sign and every aspect points the wrong way.

The weighting — double on the cells adjacent to the centre — makes Horn’s method a smoothed gradient estimate. On a noisy global DEM that smoothing is exactly what keeps a single bad cell from producing a spike of slope around it. The simpler central-difference estimate from np.gradient uses only the four direct neighbours and is noticeably noisier, which is why the version in adding terrain derivatives as predictors smooths the DEM first.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Trigonometry and array handling
scipy >=1.11 ndimage.convolve for the 3×3 kernels
rasterio >=1.3.0 Reading the DEM and writing results 
pip install "numpy>=1.23" "scipy>=1.11" "rasterio>=1.3.0"

Complete Working Example

import numpy as np
import rasterio
from scipy import ndimage

KX = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype="float64") / 8.0
KY = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]], dtype="float64") / 8.0


def slope_aspect(dem_path: str, *, flat_deg: float = 0.5):
    with rasterio.open(dem_path) as src:
        if not src.crs.is_projected:
            raise ValueError("DEM must be in a projected CRS in metres")
        dem = src.read(1, masked=True).astype("float64")
        res_x, res_y = abs(src.transform.a), abs(src.transform.e)
        profile = src.profile

    z = dem.filled(np.nan)
    valid = np.isfinite(z)
    z_filled = np.where(valid, z, np.nanmean(z))           # avoid NaN spreading via the kernel

    dzdx = ndimage.convolve(z_filled, KX, mode="nearest") / res_x
    dzdy = ndimage.convolve(z_filled, KY, mode="nearest") / res_y

    slope = np.degrees(np.arctan(np.hypot(dzdx, dzdy))).astype("float32")
    aspect = ((np.degrees(np.arctan2(-dzdx, -dzdy)) + 360) % 360).astype("float32")
    aspect[slope < flat_deg] = np.nan                       # flat: no direction

    # A cell whose 3x3 neighbourhood touches nodata is not trustworthy
    near_void = ndimage.binary_dilation(~valid, iterations=1)
    slope[near_void] = np.nan
    aspect[near_void] = np.nan
    return slope, aspect, profile

Two protective measures are worth noting. Filling voids with the mean before convolving stops NaN spreading through the kernel, and then masking every cell adjacent to a void afterwards removes the gradients that the temporary fill would otherwise have fabricated. The result is honest: slope where there is evidence, NaN where there is not. The proper treatment of voids is in filling voids in DEM tiles.

Aspect Conventions

Geographic and mathematical angles Geographic aspect is measured clockwise from north, so east is 90 degrees and south 180. The mathematical convention measures anticlockwise from east. A slope facing east is 90 degrees geographically and 0 degrees mathematically. Mixing the conventions rotates and mirrors every aspect in the raster. geographic: clockwise from north mathematical: anticlockwise from east 0° N 90° E 180° S 270° W 90° 0° 270° 180° Every GIS and every illumination formula expects the left-hand convention.

Geographic aspect, clockwise from north, is what GIS software, hillshading and illumination-correction formulas expect. np.arctan2 returns the mathematical angle, anticlockwise from east, so the arguments must be arranged to produce the geographic angle directly — which is what arctan2(-dzdx, -dzdy) does for a downslope direction. Reversing either sign silently rotates or mirrors every value, and the output still contains plausible numbers between 0 and 360.

Processing Large DEMs in Windows

A national DEM at 30 m is tens of thousands of pixels on a side, and the convolution approach above holds it all in memory. For large areas, process in windows with a one-pixel halo so the 3×3 kernel has the neighbours it needs at every window edge, then write only the window’s interior. Without the halo, every window boundary gets a row and column of edge-handled gradients that show as a faint grid in the slope layer.

The pattern is identical to the halo approach used for model inference — read one pixel more on each side than you write — and it parallelises across windows with no coordination, since each window’s output depends only on its own read. For very large extents, running the windows as independent tasks, as described in scaling raster processing with Dask, makes a continental slope layer a routine job.

Verification

Test on a surface whose answer you know A synthetic plane rising to the north at a known gradient must produce a uniform slope equal to that gradient's angle and a uniform aspect of 180 degrees, facing south and downhill. A plane rising to the east must produce aspect of 270. Any other answer means a sign or axis is wrong. Synthetic planes with known answers rises to the north, 10% slope 5.71°, aspect 180° faces south, downhill rises to the east, 10% slope 5.71°, aspect 270° faces west, downhill Two tiny arrays test every sign in the implementation. 
import numpy as np
from scipy import ndimage

res = 30.0
rows, cols = np.mgrid[0:50, 0:50]
north_plane = (49 - rows) * res * 0.10                     # row 0 is north: rises northward
dzdx = ndimage.convolve(north_plane, KX, mode="nearest") / res
dzdy = ndimage.convolve(north_plane, KY, mode="nearest") / res
s = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))[5:-5, 5:-5]
a = ((np.degrees(np.arctan2(-dzdx, -dzdy)) + 360) % 360)[5:-5, 5:-5]
assert np.allclose(s, np.degrees(np.arctan(0.10)), atol=0.01)
assert np.allclose(a, 180.0, atol=0.01), a.mean()

A synthetic plane is the definitive test because the right answer is known analytically. Keeping this check in the test suite means any later change to the kernel, the signs or the aspect formula is caught immediately rather than discovered in a hillshade that looks faintly wrong.

Common Errors

Aspect is mirrored east to west

The sign of the east gradient is inverted. Test on a plane rising to the east.

North-facing slopes are reported as south-facing

The north gradient ignores that row zero is northernmost. Flip the kernel’s vertical sign.

A ring of steep slope around every void

NaN or fill values entered the kernel. Fill temporarily, then mask cells adjacent to voids.

Flat areas show a patchwork of aspects

Aspect was computed for zero-gradient cells. Flag cells below a small slope threshold as undefined.

Frequently Asked Questions

Q: Horn’s method or simple central differences? Horn’s method, for most purposes. It uses all eight neighbours with a weighting that suppresses noise, which matters on the noisy surfaces of global DEMs. Simple central differences use only four neighbours and are more sensitive to single-cell errors.

Q: What is the aspect of a perfectly flat cell? Undefined. With zero gradient there is no downslope direction, and any number the formula returns is arbitrary. Flag flat cells explicitly, commonly with a sentinel value or NaN, rather than letting them default to zero degrees, which means north.

Q: Should slope be in degrees or per cent? Either, as long as it is stated. Degrees are conventional in analysis and modelling; per cent rise is common in engineering and agriculture. They are related by the tangent: 45 degrees is 100 per cent.

Q: Does GDAL’s DEM tool give the same answer? Yes, when configured with Horn’s algorithm, which is its default, and the same edge handling. It is a convenient cross-check: run it on the same DEM and compare, and any difference beyond floating-point noise points at a sign or scale error.

Q: How should aspect be used as a model feature? As its sine and cosine rather than as degrees, because aspect is circular — 359 and 1 degrees face the same way. The encoding is set out in adding terrain derivatives as predictors.