How to Calculate Slope and Aspect from a DEM in Python

On this page

Problem statement

There is no rasterio.slope(). Slope and aspect are derived products you compute yourself, or shell out to GDAL for, and both routes have a trap that produces plausible nonsense:

import numpy as np
import rasterio

with rasterio.open("dem.tif") as src:
    dem = src.read(1)
    dy, dx = np.gradient(dem, src.res[1], src.res[0])

slope = np.degrees(np.arctan(np.hypot(dx, dy)))
print(f"mean {slope.mean():.1f}°  max {slope.max():.1f}°")

mean 88.6°  max 90.0°

Snowdonia is not, on average, 88.8° steep. src.res is in degrees because the DEM is in EPSG:4326, so the gradient is 67,000 to 110,000 times too large and arctan saturates at vertical.

Three things have to be right: the horizontal cell size in ground units, NoData masked before the window arithmetic, and the aspect convention. This guide does all three, and shows the GDAL route for when you would rather not.

Quick answer

import math

import numpy as np
import rasterio

with rasterio.open("dem.tif") as src:
    dem = src.read(1).astype("float64")
    if src.nodata is not None:
        dem = np.where(dem == src.nodata, np.nan, dem)

    if src.crs.is_geographic:
        lat = (src.bounds.bottom + src.bounds.top) / 2
        cell_x = abs(src.transform.a) * 111_320 * math.cos(math.radians(lat))
        cell_y = abs(src.transform.e) * 110_574
    else:
        cell_x, cell_y = abs(src.transform.a), abs(src.transform.e)

p = np.pad(dem, 1, mode="edge")
a, b, c = p[:-2, :-2], p[:-2, 1:-1], p[:-2, 2:]
d,    f = p[1:-1, :-2],               p[1:-1, 2:]
g, h, i = p[2:, :-2],  p[2:, 1:-1],  p[2:, 2:]

dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cell_x)      # Horn
dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cell_y)

gradient = np.hypot(dz_dx, dz_dy)
slope = np.degrees(np.arctan(gradient))
aspect = np.where(gradient < 1e-8, np.nan, np.degrees(np.arctan2(dz_dy, -dz_dx)) % 360)

print(f"cell {cell_x:.1f} x {cell_y:.1f} m")
print(f"slope mean {np.nanmean(slope):.2f}°  max {np.nanmax(slope):.2f}°")

cell 27.9 x 30.7 m
slope mean 21.31°  max 65.70°

Step two is where the catastrophic errors live, and it is the step with no code in most tutorials.

Mask, ground spacing, Horn's window, sanity checks

1. Read and mask NoData first

with rasterio.open("dem.tif") as src:
    dem = src.read(1).astype("float64")
    if src.nodata is not None:
        dem = np.where(dem == src.nodata, np.nan, dem)
dem = np.where(dem < -100, np.nan, dem)        # undeclared sentinels: -32768, -9999
print(f"{np.isnan(dem).sum()} NoData cells")

0 NoData cells

Do this before any window arithmetic. This Copernicus window has none, but a single -32768 — SRTM's void value — puts a slope of 89° or more on the eight cells around it, so a NoData region gets a one-cell halo of fake cliffs around its whole boundary.

.astype("float64") matters too: integer DEMs are common, and integer division truncates the gradient.

2. Get the cell size in ground units

with rasterio.open("dem.tif") as src:
    if src.crs.is_geographic:
        lat = (src.bounds.bottom + src.bounds.top) / 2
        cell_x = abs(src.transform.a) * 111_320 * math.cos(math.radians(lat))
        cell_y = abs(src.transform.e) * 110_574
    else:
        cell_x, cell_y = abs(src.transform.a), abs(src.transform.e)
print(f"{cell_x:.1f} x {cell_y:.1f} m")

27.9 x 30.7 m

Two things worth noting. transform.e is negative for a north-up raster, so take the absolute value. And do not assume the two are equal — Copernicus DEM at 53°N is 27.9 × 30.7 m, and using cell_x for both introduces a 4.5% error in the mean slope.

The alternative, when accuracy matters more than convenience, is to reproject the DEM to a metric CRS first — but that resamples the elevations, which introduces its own smoothing. Deriving the cell size is usually the better trade.

3. Apply Horn's 3×3 window

p = np.pad(dem, 1, mode="edge")
a, b, c = p[:-2, :-2], p[:-2, 1:-1], p[:-2, 2:]
d,    f = p[1:-1, :-2],               p[1:-1, 2:]
g, h, i = p[2:, :-2],  p[2:, 1:-1],  p[2:, 2:]

dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cell_x)
dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cell_y)

The nine slices are the 3×3 neighbourhood of every cell at once — fully vectorised, no Python loop. On a 252 × 240 grid this runs in milliseconds; on a 10,000² DEM it is still seconds rather than hours.

Horn weights the direct neighbours twice and the diagonals once. It is what GDAL, QGIS and ArcGIS use, so your numbers will match theirs.

4. Convert to slope and aspect

gradient = np.hypot(dz_dx, dz_dy)

slope_deg = np.degrees(np.arctan(gradient))
slope_pct = gradient * 100
aspect = np.where(gradient < 1e-8, np.nan, np.degrees(np.arctan2(dz_dy, -dz_dx)) % 360)

The -dz_dx is not a typo: aspect points downslope, the direction water flows. And the where guard matters — on flat ground arctan2(0, 0) returns 0, which reads as due north and produces a large phantom north-facing region.

5. Sanity-check before writing

if np.nanmax(slope_deg) > 89:
    raise ValueError(f"max slope {np.nanmax(slope_deg):.1f}° — check cell units and NoData")
print(f"over 45°: {np.nanmean(slope_deg > 45):.1%}")

over 45°: 2.8%

Real terrain over a 30 m cell almost never exceeds 70°. A maximum near 90° is a bug every time, and this one line catches both the units error and NoData contamination.

The same cell in two unit systems. The left-hand numbers produce a gradient 67,000 to 110,000 times too large.

A reusable terrain function and the GDAL route

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, 0.1° × 0.07° around the summit, read from the public bucket and saved twice under the two names this page uses. The window crosses the 4°W tile boundary, so two tiles are merged.

import shutil
import rasterio
from rasterio.merge import merge

BBOX = (-4.10, 53.03, -4.00, 53.10)
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
       "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
tiles = [COP.format("W005"), COP.format("W004")]
mosaic, transform = merge(tiles, bounds=BBOX)
with rasterio.open(tiles[0]) as src:
    profile = src.profile | {"driver": "GTiff", "width": mosaic.shape[2], "height": mosaic.shape[1],
                             "transform": transform, "compress": "deflate"}
with rasterio.open("dem.tif", "w", **profile) as dst:
    dst.write(mosaic)
shutil.copy("dem.tif", "snowdonia_glo30.tif")
print(f"dem.tif {mosaic.shape[2]} × {mosaic.shape[1]}, {mosaic.min():.1f} to {mosaic.max():.1f} m")

dem.tif 240 × 252, 54.0 to 1074.5 m

Example 1 — a complete, reusable terrain function

import math

import numpy as np
import rasterio


def cell_size_metres(src):
    """Ground cell size in metres for either a projected or a geographic raster."""
    if not src.crs.is_geographic:
        return abs(src.transform.a), abs(src.transform.e)
    lat = (src.bounds.bottom + src.bounds.top) / 2
    return (abs(src.transform.a) * 111_320 * math.cos(math.radians(lat)),
            abs(src.transform.e) * 110_574)


def slope_aspect(path, *, method="horn", nodata_below=-100, max_plausible=80.0):
    """Slope (degrees) and aspect (degrees from north) from a DEM."""
    with rasterio.open(path) as src:
        dem = src.read(1).astype("float64")
        if src.nodata is not None:
            dem = np.where(dem == src.nodata, np.nan, dem)
        cell_x, cell_y = cell_size_metres(src)
        profile = src.profile.copy()
        geographic = src.crs.is_geographic

    invalid = ~np.isfinite(dem) | (dem < nodata_below)
    if invalid.any():
        print(f"  masked {invalid.sum():,} invalid cell(s)")
    working = np.where(invalid, np.nanmean(dem[~invalid]), dem)

    p = np.pad(working, 1, mode="edge")
    a, b, c = p[:-2, :-2], p[:-2, 1:-1], p[:-2, 2:]
    d,    f = p[1:-1, :-2],               p[1:-1, 2:]
    g, h, i = p[2:, :-2],  p[2:, 1:-1],  p[2:, 2:]

    if method == "horn":
        dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cell_x)
        dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cell_y)
    elif method == "zevenbergen":
        dz_dx = (f - d) / (2 * cell_x)
        dz_dy = (h - b) / (2 * cell_y)
    else:
        raise ValueError(f"unknown method {method!r}")

    gradient = np.hypot(dz_dx, dz_dy)
    slope = np.degrees(np.arctan(gradient))
    aspect = np.where(gradient < 1e-8, np.nan,
                      np.degrees(np.arctan2(dz_dy, -dz_dx)) % 360)

    # grow the invalid mask by one cell — the window contaminates its neighbours
    from scipy.ndimage import binary_dilation
    halo = binary_dilation(invalid, np.ones((3, 3), bool))
    slope = np.where(halo, np.nan, slope)
    aspect = np.where(halo, np.nan, aspect)

    print(f"  {'geographic' if geographic else 'projected'} · cell {cell_x:.1f} x {cell_y:.1f} m · {method}")
    print(f"  slope mean {np.nanmean(slope):.2f}°  median {np.nanmedian(slope):.2f}°  "
          f"max {np.nanmax(slope):.2f}°  over 45°: {np.nanmean(slope > 45):.1%}")

    if np.nanmax(slope) > max_plausible:
        raise ValueError(
            f"max slope {np.nanmax(slope):.1f}° exceeds {max_plausible}° — "
            f"cell size or NoData is wrong"
        )
    return slope, aspect, profile, (cell_x, cell_y)


slope, aspect, profile, cell = slope_aspect("snowdonia_glo30.tif")

  geographic · cell 27.9 x 30.7 m · horn
  slope mean 21.31°  median 20.32°  max 65.70°  over 45°: 2.8%

Two details worth keeping. The invalid mask is dilated by one cell before being applied to the output, because a 3×3 window spreads contamination outward — masking only the original cells leaves a ring of fake cliffs. And the max_plausible check raises rather than warns, because a slope raster with 90° values will silently poison every downstream product.

Example 2 — writing outputs that stay interpretable

def write_terrain(slope, aspect, profile, prefix, *, cell, method):
    profile.update(dtype="float32", count=1, nodata=np.nan,
                   compress="deflate", tiled=True, blockxsize=256, blockysize=256)

    outputs = [
        ("slope_deg", slope, "degrees"),
        ("slope_pct", np.tan(np.radians(slope)) * 100, "percent"),
        ("aspect", aspect, "degrees_clockwise_from_north"),
    ]
    for name, data, unit in outputs:
        path = f"{prefix}_{name}.tif"
        with rasterio.open(path, "w", **profile) as dst:
            dst.write(data.astype("float32"), 1)
            dst.update_tags(quantity=name, unit=unit, method=method,
                            cell_x_m=f"{cell[0]:.2f}", cell_y_m=f"{cell[1]:.2f}")
        print(f"  {path:28} {unit}")


write_terrain(slope, aspect, profile, "snowdonia", cell=cell, method="horn")

  snowdonia_slope_deg.tif      degrees
  snowdonia_slope_pct.tif      percent
  snowdonia_aspect.tif         degrees_clockwise_from_north

Writing both slope units removes an entire class of downstream confusion: a value of 45 means very different things in the two files, and only the tags distinguish them.

Every bar is a line printed by Example 3 below or by Example 1.

Example 3 — the GDAL route, and why its scale factor cannot fix a geographic DEM

gdaldem slope has a -s option described as "ratio of vertical units to horizontal", with 111120 offered as the value for a lat/long DEM. On this DEM none of the obvious choices is right:

import subprocess

import rasterio


def gdaldem_slope(src_path, dst_path, *, scale=1.0):
    subprocess.run(["gdaldem", "slope", src_path, dst_path,
                    "-s", str(scale), "-q",
                    "-co", "COMPRESS=DEFLATE"], check=True)
    with rasterio.open(dst_path) as src:
        data = src.read(1, masked=True)
    implied_x = abs(src.transform.a) * scale
    implied_y = abs(src.transform.e) * scale
    print(f"  -s {scale:<8.0f} implies cell {implied_x:5.1f} x {implied_y:5.1f} m  ->  "
          f"mean {data.mean():6.2f}°  max {data.max():6.2f}°")
    return data


for scale in (1, 111_120, 111_120 * math.cos(math.radians(53.065))):
    gdaldem_slope("snowdonia_glo30.tif", f"gdal_slope_{int(scale)}.tif", scale=scale)

print(f"  true ground cell is {cell[0]:.1f} x {cell[1]:.1f} m  ->  numpy gives mean {np.nanmean(slope):.2f}°  max {np.nanmax(slope):.2f}°")

  -s 1        implies cell   0.0 x   0.0 m  ->  mean  88.78°  max  90.00°
  -s 111120   implies cell  46.3 x  30.9 m  ->  mean  18.08°  max  64.66°
  -s 66773    implies cell  27.8 x  18.5 m  ->  mean  27.35°  max  74.12°
  true ground cell is 27.9 x 30.7 m  ->  numpy gives mean 21.31°  max 65.70°

Read the middle column. -s multiplies both axes by the same number, so:

111120 makes the y cell 30.9 m — correct — and the x cell 46.3 m, 66% too large. Slopes come out 15% too shallow.
111120 × cos(lat), the value people reach for to "fix the longitude", makes x correct at 27.8 m and shrinks y to 18.5 m. Slopes come out 28% too steep.

Neither is right, because a Copernicus DEM's ground cell is 27.9 × 30.7 m and no single multiplier produces both from 0.000417 and 0.000278 degrees.

There are exactly two correct options:

# a) reproject to a metric CRS first, then -s 1 is genuinely correct

for extra, name in [([], "no_nodata"), (["-dstnodata", "-9999"], "nodata")]:
    subprocess.run(["gdalwarp", "-q", "-overwrite", "-t_srs", "EPSG:27700", "-r", "bilinear", *extra,
                    "snowdonia_glo30.tif", f"snowdonia_bng_{name}.tif"], check=True)
    gdaldem_slope(f"snowdonia_bng_{name}.tif", f"slope_bng_{name}.tif", scale=1.0)

# b) compute it yourself with per-axis cell sizes, as in Example 1

  -s 1        implies cell  29.5 x  29.5 m  ->  mean  22.19°  max  86.80°
  -s 1        implies cell  29.5 x  29.5 m  ->  mean  21.19°  max  64.06°

The first reprojection is a warning in itself. The Copernicus DEM declares no NoData, so gdalwarp fills the corners of the rotated grid with 0, and the drop from mountain to "sea level" at every edge reads as an 86.8° cliff. Declared with -dstnodata, the same surface gives a mean of 21.19° and a maximum of 64.06° — close to the 21.31° and 65.70° computed on the original grid, and slightly lower, because reprojecting resamples the elevations and smooths the surface. Computing it yourself avoids that. Both beat any single -s value on a DEM whose ground spacings differ.

np.gradient, dilated masks and computing before reprojecting

Why `np.gradient` is not quite the right tool

np.gradient computes central differences on the four direct neighbours — mathematically it is close to Zevenbergen–Thorne, not Horn. It is a perfectly reasonable estimator and it will not match GDAL, QGIS or ArcGIS output.

It also has an argument-order trap:

dy, dx = np.gradient(dem, cell_y, cell_x)      # rows first, then columns

Rows are the y axis, columns are x, and np.gradient returns them in array-axis order. Swapping them transposes your aspect map by 90° while leaving slope unchanged — so slope looks fine and aspect is wrong, which is very hard to spot.

The explicit nine-slice version is more code and has no ambiguity about which axis is which.

Why the invalid mask must be dilated

Every cell's slope is computed from a 3×3 window. If any cell in that window was invalid, the centre cell's value is contaminated — even though the centre cell itself had a perfectly good elevation.

So the output mask should be the input mask grown by one cell in every direction. Skipping the dilation leaves a one-pixel ring of extreme values around every NoData region, which then survives into contours, hillshade, thresholds and any zonal statistic. It is the same class of problem as NoData in a zonal sum, one step removed.

Why replacing NaN with the mean is safe here

Example 1 substitutes the mean elevation into invalid cells before the window arithmetic, then masks the output. That looks like it should corrupt things, and it does not, because every cell whose window touched an invalid value is masked afterwards by the dilated halo.

The substitution exists purely to keep the arithmetic finite — NaN propagates through + and *, so a single NaN would poison a much larger area than the dilation covers, and some operations would produce warnings. Filling then masking is faster and clearer than a NaN-aware convolution.

Every window touching a bad cell is contaminated. The output mask must be one cell larger than the input mask.

Why to compute rather than reproject

The obvious alternative to deriving the cell size is to reproject the DEM into a metric CRS, where the cell size is already metres. It works, and it costs something: reprojection resamples the elevation values, which smooths the surface and slightly reduces the slopes you then compute.

For a small area the difference is minor and reprojecting is simpler. For a large area — where a single projected CRS distorts across the extent — deriving the cell size per latitude band is more accurate and avoids touching the elevations at all.

Either is defensible. Doing neither is not.

Negative transforms, the -s factor and integer DEMs

transform.e is negative for a north-up raster. Use abs().
gdaldem's -s applies one factor to both axes, so it cannot represent a DEM whose x and y ground spacings differ. Reproject first, or compute the gradient yourself. Without -s at all, the output saturates at 90°.
np.gradient takes spacings in array-axis order — rows (y) first. Getting it wrong rotates aspect and leaves slope unchanged.
Integer DEMs truncate. Cast to float64 before any arithmetic.
The outer ring is artificially flat with mode="edge" padding. Compute on a padded extent and crop if the boundary matters.
Aspect on flat ground is undefined. Guard on the gradient magnitude or you get a phantom north-facing plateau.
Slope from a DSM has a cliff at every forest edge. For terrain questions use a DTM — see DEM, DSM and DTM.
richdem and whitebox provide many more terrain indices (curvature, TPI, TWI) if you need them; the two here are the ones everything else builds on.

Internal links

Slope and aspect explained — what the numbers mean and why the choices matter
Slope values are wrong or absurdly steep — diagnosing a broken slope raster
Digital elevation models explained — the DEM this all starts from
How to make a hillshade from a DEM in Python — slope and aspect combined into an image
The raster data model explained — the transform the cell size comes from
How to calculate zonal statistics in Python — summarising slope by region
Raster resampling explained — the cost of reprojecting a DEM first
How to generate contour lines from a DEM in Python — the other classic derivative
Landslide susceptibility explained — slope as a hazard factor, and what it cannot say alone
Terrain shading and the horizon explained — what slope and aspect do to solar yield

FAQ

Why is my slope 90° everywhere?

The horizontal cell size is in degrees. Convert to metres using the latitude before computing the gradient.

Can I just use `np.gradient`?

You can, but it implements central differences rather than Horn's method, so it will not match GDAL or QGIS. It also takes spacings in row-then-column order, which is easy to reverse.

How do I get the cell size in metres from a geographic DEM?

Multiply the longitude spacing by 111,320 × cos(latitude) and the latitude spacing by 110,574. Do not assume the two results are equal.

Should I reproject the DEM first instead?

Either works. Reprojecting resamples the elevations and smooths the surface slightly; deriving the cell size leaves them untouched. For large areas, deriving is more accurate.

Why is there a ring of steep cells around my NoData?

Every 3×3 window touching a bad cell is contaminated. Dilate the invalid mask by one cell before applying it to the output.

What does `gdaldem slope -s` do?

It scales the horizontal units to match the vertical unit — but with a single factor for both axes. On a Copernicus DEM, whose x and y ground spacings differ, no value is correct: 111120 understates slope by 15% and 111120·cos(lat) overstates it by 28%. Reproject to a metric CRS first, or compute the gradient yourself.

How do I convert degrees to percent?

tan(radians(slope)) * 100. 45° is 100%, and percent has no upper bound.