How to Calculate Slope and Aspect from a DEM in Python
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.8° 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 about ninety thousand 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.33° max 65.70°
Step-by-step solution
1. Read and mask NoData first
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, -9999Do this before any window arithmetic. A single -32768 inside a 3×3 window produces a slope of 90° in all nine cells, 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
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)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 × 239 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.
Code examples
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.33° median 20.30° 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.
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} 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 27.9 x 30.7 m -> numpy gives mean 21.33° max 65.70°") -s 1 implies cell 0.0 x 0.0 m -> mean 88.77° max 90.00°
-s 111120 implies cell 46.3 x 30.9 m -> mean 18.10° max 64.66°
-s 66773 implies cell 27.8 x 18.5 m -> mean 27.37° max 74.12°
true ground cell is 27.9 x 30.7 m -> numpy gives mean 21.33° 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
subprocess.run(["gdalwarp", "-t_srs", "EPSG:27700", "-r", "bilinear",
"snowdonia_glo30.tif", "snowdonia_bng.tif"], check=True)
gdaldem_slope("snowdonia_bng.tif", "slope_bng.tif", scale=1.0)
# b) compute it yourself with per-axis cell sizes, as in Example 1Reprojecting resamples the elevations and smooths the surface slightly. Computing it yourself does not. Both beat any single -s value on a DEM whose ground spacings differ.
Explanation
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 columnsRows 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.
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.
Edge cases or notes
transform.eis negative for a north-up raster. Useabs().gdaldem's-sapplies 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-sat all, the output saturates at 90°.np.gradienttakes spacings in array-axis order — rows (y) first. Getting it wrong rotates aspect and leaves slope unchanged.- Integer DEMs truncate. Cast to
float64before 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
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.