🧭 Concept Terrain & elevation Raster & rasterio Terrain & Elevation · step 3 of 12 Updated

Slope and Aspect Explained: How They Are Computed from a Grid

Problem statement

Slope looks like it should be simple: how steep is the ground. But a DEM is a grid of numbers, and "steepness at a point" is not something a grid stores — it has to be estimated from the neighbours.

That estimate involves at least four decisions, and every one of them changes the answer:

  • Which neighbours — 4 or 8? Weighted how?
  • What horizontal distance do the neighbours sit at? On a geographic DEM this is not the cell size.
  • What units — degrees, percent, or a ratio?
  • What happens at the edges and at NoData?

Get the second one wrong and the numbers are not slightly off, they are nonsense:

                                     mean      max   over 45°
degrees used as the cell size      88.792°  90.000°     98.7%
metres                             21.326°  65.698°      2.8%

Same DEM. The first row says 98.7% of Snowdonia is steeper than 45°, which would make it unwalkable. Nothing raised.

Quick answer

Slope is the magnitude of the elevation gradient; aspect is its direction.

import numpy as np


def horn_gradient(dem, cell_x, cell_y):
    """Horn's method — the 3x3 weighted estimate used by GDAL, ArcGIS and QGIS."""
    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)
    return dz_dx, dz_dy


dz_dx, dz_dy = horn_gradient(dem, cell_x=27.9, cell_y=30.7)     # metres, not degrees

slope_deg = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
aspect_deg = np.degrees(np.arctan2(dz_dy, -dz_dx)) % 360        # 0 = north, clockwise

print(f"slope  mean {slope_deg.mean():.2f}°  max {slope_deg.max():.2f}°")
slope  mean 21.33°  max 65.70°
Quantity Formula Range
gradient magnitude hypot(dz_dx, dz_dy) 0 to ∞ (rise over run)
slope in degrees degrees(arctan(gradient)) 0 to 90
slope in percent gradient * 100 0 to ∞
aspect degrees(arctan2(dz_dy, -dz_dx)) % 360 0 to 360, 0 = north
A three by three window of elevation cells with Horn's weights, showing how the east-west and north-south gradients are formed.
Horn's method weights the direct neighbours twice and the diagonals once, in each direction.

Step-by-step solution

1. Get the cell size in ground units

This is the step that produces the catastrophic errors. A DEM in EPSG:4326 has a cell size in degrees:

print(src.res)
(0.00041666666666666664, 0.0002777777777777778)

Passing those numbers as the horizontal distance means dividing a height difference in metres by a distance of 0.0004 — a gradient about 90,000 times too large, which arctan then saturates at 90°.

Convert first:

import math

lat = 53.065
cell_x = src.res[0] * 111_320 * math.cos(math.radians(lat))
cell_y = src.res[1] * 110_574
print(f"{cell_x:.1f} m x {cell_y:.1f} m")
27.9 m x 30.7 m

Note they differ. Assuming square cells and using cell_x for both gives a mean slope of 22.28° against the correct 21.33° — a 4.5% error, which is quiet enough to survive review.

2. Choose the gradient estimator

Two are in common use:

Horn (1981) — the 3×3 weighted estimate above. Used by GDAL, ArcGIS and QGIS. It weights the four direct neighbours twice and the four diagonals once, which smooths noise a little.

Zevenbergen–Thorne (1987) — fits a partial quartic and uses only the four direct neighbours:

def zevenbergen_thorne(dem, cell_x, cell_y):
    p = np.pad(dem, 1, mode="edge")
    dz_dx = (p[1:-1, 2:] - p[1:-1, :-2]) / (2 * cell_x)
    dz_dy = (p[2:, 1:-1] - p[:-2, 1:-1]) / (2 * cell_y)
    return dz_dx, dz_dy

Horn is smoother and is the default nearly everywhere; Zevenbergen–Thorne preserves sharp breaks better and is noisier. On the DEM used here the two differ by 0.15° in the mean — real, but far smaller than a units mistake. Pick one, and record which.

3. Choose the units and mean it 

gradient = np.hypot(dz_dx, dz_dy)

slope_deg = np.degrees(np.arctan(gradient))
slope_pct = gradient * 100
slope_ratio = gradient

Degrees are bounded at 90 and intuitive. Percent is unbounded — a vertical cliff is infinite percent — and is what road and rail engineering uses. A "10% slope" is 5.7°, and people confuse the two constantly.

4. Understand aspect's conventions

Aspect is a compass bearing: 0° is north, increasing clockwise. Two details trip people up.

The sign of dz_dx. Aspect points downslope — the direction water flows — so the x-gradient is negated:

aspect = np.degrees(np.arctan2(dz_dy, -dz_dx)) % 360

Flat cells have no aspect. Where the gradient is zero, arctan2(0, 0) returns 0, which reads as "north". Flag it:

aspect = np.where(gradient < 1e-8, np.nan, aspect)

Without that, a flat plateau appears as a large north-facing area — the commonest artefact in an aspect map.

5. Handle edges and NoData 

p = np.pad(dem, 1, mode="edge")

mode="edge" replicates the boundary row, giving zero gradient across it — so the outer ring reads as flatter than reality. Alternatives are to compute on a padded extent and crop, or to set the outer ring to NaN and be explicit about it.

NoData is worse. A single -32768 inside the window produces a gradient of about 1,000,000 and a slope of 90° in all eight surrounding cells. Mask it before computing anything.

The same gradient expressed as degrees, percent and a ratio, showing 45 degrees equals 100 percent.
45° is 100%, not 50%. The two scales diverge fast, and percent has no upper bound.

Code examples

Example 1 — slope and aspect with the cell size derived, not assumed 

import math

import numpy as np
import rasterio


def cell_size_metres(src):
    """Ground cell size in metres, correct for geographic and projected CRSs."""
    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 terrain(path, *, method="horn", nodata_below=-100):
    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)
        dem = np.where(dem < nodata_below, np.nan, dem)
        cell_x, cell_y = cell_size_metres(src)
        profile = src.profile.copy()

    print(f"  cell {cell_x:.1f} m x {cell_y:.1f} m ({'geographic' if src.crs.is_geographic else 'projected'})")

    filled = np.where(np.isnan(dem), np.nanmean(dem), dem)   # keep the window arithmetic finite
    p = np.pad(filled, 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)
    else:                                                     # zevenbergen-thorne
        dz_dx = (f - d) / (2 * cell_x)
        dz_dy = (h - b) / (2 * 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)

    invalid = np.isnan(dem)
    slope = np.where(invalid, np.nan, slope)
    aspect = np.where(invalid, np.nan, aspect)

    if np.nanmax(slope) > 89:
        print(f"  WARNING: max slope {np.nanmax(slope):.1f}° — check cell units and NoData")
    print(f"  slope: mean {np.nanmean(slope):.2f}°  max {np.nanmax(slope):.2f}°  "
          f"over 45°: {np.nanmean(slope > 45):.1%}")
    return slope, aspect, profile


slope, aspect, profile = terrain("snowdonia_glo30.tif")
  cell 27.9 m x 30.7 m (geographic)
  slope: mean 21.33°  max 65.70°  over 45°: 2.8%

The warning threshold catches the unit error every time: a real landscape has almost nothing above 70°, so a max near 90° means the horizontal distance is wrong.

Example 2 — what the four decisions actually cost 

def compare_methods(path):
    import pandas as pd

    with rasterio.open(path) as src:
        dem = src.read(1).astype("float64")
        res_x, res_y = abs(src.transform.a), abs(src.transform.e)
        cell_x, cell_y = cell_size_metres(src)

    def slope_from(xr, yr, method="horn"):
        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:]
        if method == "horn":
            gx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * xr)
            gy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * yr)
        else:
            gx, gy = (f - d) / (2 * xr), (h - b) / (2 * yr)
        return np.degrees(np.arctan(np.hypot(gx, gy)))

    rows = [
        ("degrees as cell size", slope_from(res_x, res_y)),
        ("square cells (x for both)", slope_from(cell_x, cell_x)),
        ("correct metres, Horn", slope_from(cell_x, cell_y)),
        ("correct metres, Z-T", slope_from(cell_x, cell_y, "zt")),
    ]
    frame = pd.DataFrame([
        {"variant": name, "mean": round(s.mean(), 3), "max": round(s.max(), 3),
         "over_45": f"{(s > 45).mean():.1%}"}
        for name, s in rows
    ])
    print(frame.to_string(index=False))


compare_methods("snowdonia_glo30.tif")
                  variant   mean    max over_45
     degrees as cell size 88.792 90.000   98.7%
square cells (x for both) 22.284 67.483    3.9%
     correct metres, Horn 21.326 65.698    2.8%
      correct metres, Z-T 21.479 66.643    3.1%

Row one is unmistakably broken. Row two is the dangerous one — 22.28° against 21.33° is a 4.5% error in the mean and a 39% error in "how much land is steeper than 45°" (3.9% against 2.8%), from an assumption nobody writes down.

Rows three and four differ by 0.15° in the mean, which is the honest size of the Horn versus Zevenbergen–Thorne choice on a clean DEM: real, small, and worth recording rather than agonising over.

Example 3 — writing the outputs so they stay interpretable 

def write_terrain(slope, aspect, profile, prefix, *, cell_x, cell_y, method):
    profile.update(dtype="float32", count=1, nodata=np.nan, compress="deflate")

    for name, data, unit in [("slope", slope, "degrees"), ("aspect", aspect, "degrees_from_north")]:
        with rasterio.open(f"{prefix}_{name}.tif", "w", **profile) as dst:
            dst.write(data.astype("float32"), 1)
            dst.update_tags(
                quantity=name, unit=unit, method=method,
                cell_x_m=f"{cell_x:.2f}", cell_y_m=f"{cell_y:.2f}",
                note="aspect is NaN where the surface is flat",
            )
        print(f"  wrote {prefix}_{name}.tif ({unit}, {method})")


write_terrain(slope, aspect, profile, "snowdonia",
              cell_x=27.9, cell_y=30.7, method="horn")
  wrote snowdonia_slope.tif (degrees, horn)
  wrote snowdonia_aspect.tif (degrees_from_north, horn)

A slope raster without its units recorded is ambiguous between degrees, percent and a ratio — three interpretations that differ by orders of magnitude. Four tags cost nothing.

Explanation

Why slope is scale-dependent

Slope is not a property of a point. It is a property of a point at a given resolution, because it is estimated over a 3×3 window whose size is the cell size.

Compute slope from a 1 m LiDAR DTM and you get individual boulders and kerbs. Compute it from a 30 m DEM over the same hillside and you get the average of a 90 m patch. Both are correct; they measure different things.

Resampling a 30 m DEM to 5 m before computing slope does not recover detail. It produces a smoother slope map with more cells, which reads as higher quality and contains no more information. Always report the source resolution alongside a slope statistic.

Why aspect is circular and how that breaks statistics

Aspect is an angle, so 359° and 1° are 2° apart, not 358°. Ordinary statistics do not survive this:

angles = np.array([359.0, 1.0])
print(f"arithmetic mean: {angles.mean():.1f}°")

radians = np.radians(angles)
circular = np.degrees(np.arctan2(np.sin(radians).mean(), np.cos(radians).mean())) % 360
print(f"circular mean:   {circular:.1f}°")
arithmetic mean: 180.0°
circular mean:   0.0°

180° — due south — for two cells that both face almost due north. Any mean, standard deviation or interpolation of aspect must use circular statistics, or be done on the sine and cosine components separately.

Bin into compass octants rather than averaging, wherever you can.

Two aspects at 359 and 1 degrees averaging arithmetically to 180 degrees due south, and circularly to 0 degrees due north.
Two north-facing cells averaging to due south. Aspect is circular, and ordinary statistics do not apply.

Why the horizontal distance matters so much

The gradient is a rise over a run. The rise is in metres of elevation; the run must be in metres of ground distance.

On a geographic DEM the run is expressed in degrees — about 0.0004 — while the rise is tens of metres. The ratio is roughly 90,000 times too large, arctan of anything that big is essentially 90°, and every cell saturates.

The subtler failure is the square-cell assumption. Copernicus DEM varies its longitude spacing by latitude band precisely so ground cells stay near-square, but "near" is not "exactly": at 53°N it is 27.9 × 30.7 m, a 10% difference. Carrying that into the gradient produces the 4.5% error in row two of the comparison.

Why NoData contaminates a neighbourhood, not a cell

Every terrain derivative reads a 3×3 window. One bad value inside that window corrupts the centre cell — so a single NoData pixel produces a 3×3 block of nonsense, and a NoData region produces a one-cell halo of extreme slopes around its entire boundary.

Those halos look like cliffs and survive into contours, hillshades and any threshold you apply. Mask NoData to NaN before computing, and mask the derivative afterwards, as in Example 1.

Edge cases or notes

  • Percent slope is unbounded. 45° is 100%, 60° is 173%, vertical is infinite.
  • Aspect is undefined where the surface is flat. Set it to NaN or a sentinel; arctan2(0, 0) returns 0, which reads as north.
  • Aspect cannot be averaged or interpolated arithmetically. Use circular statistics or bin into octants.
  • np.pad(mode="edge") makes the outer ring artificially flat. Compute on a padded extent and crop where accuracy at the boundary matters.
  • GDAL's gdaldem slope has a -s scale option for exactly the geographic-DEM case: -s 111120 converts degrees to metres approximately.
  • A max slope near 90° is a bug, not terrain. Real landscapes rarely exceed 70° over a 30 m cell.
  • Horn and Zevenbergen–Thorne differ by a fraction of a degree on a clean DEM, and more where the surface is noisy. Record which you used.
  • Slope from a DSM has cliffs at every forest edge. Use a DTM for terrain questions — see DEM, DSM and DTM.

FAQ

What is the difference between slope in degrees and percent?

Percent is the gradient times 100, degrees is its arctangent. 45° equals 100%, not 50%, and percent has no upper bound.

Why is my maximum slope exactly 90°?

The horizontal cell size is wrong — almost always degrees being used where metres are needed. arctan saturates and every steep cell reads as vertical.

Which method should I use, Horn or Zevenbergen–Thorne?

Horn, unless you have a reason otherwise. It is the default in GDAL, QGIS and ArcGIS, so your results will match other tools. The difference is around 0.15° in the mean on a clean DEM — much smaller than getting the cell size wrong.

Why does my aspect map have a large north-facing area?

Flat cells. arctan2(0, 0) returns 0, which reads as due north. Set aspect to NaN where the gradient is essentially zero.

Can I average aspect values?

Not arithmetically — 359° and 1° average to 180°. Use circular statistics, or bin into compass octants.

Does slope depend on the DEM resolution?

Yes, strongly. A 1 m DTM gives steeper, more detailed slopes than a 30 m DEM over the same ground. Always report the source resolution.

What happens at NoData?

One bad cell corrupts every 3×3 window containing it, producing a halo of extreme slopes. Mask NoData to NaN before computing and mask the output afterwards.