🔧 Fix Terrain & elevation CRS & projections Terrain & Elevation · step 5 of 12 Updated

Slope Values Are Wrong or Absurdly Steep

Problem statement

The slope raster came out and the numbers cannot be right:

print(f"mean {slope.mean():.2f}°  max {slope.max():.2f}°  over 45°: {(slope > 45).mean():.1%}")
mean 88.79°  max 90.00°  over 45°: 98.7%

Ninety-eight percent of the landscape steeper than 45°. Snowdonia is steep; it is not a wall.

Or the failure is quieter — everything looks plausible and is 15% off. Or there is a ring of 90° cells around every lake. Or the slope map is fine and the aspect map has an enormous north-facing plateau that does not exist.

Slope is a derived product, and every one of these is an input problem rather than a formula problem.

Quick answer

Print the four numbers that determine the answer:

import math

with rasterio.open("dem.tif") as src:
    print(f"CRS          {src.crs} (geographic: {src.crs.is_geographic})")
    print(f"res          {src.res}")
    print(f"nodata       {src.nodata}")
    dem = src.read(1)
    print(f"dtype        {dem.dtype}")
    print(f"elev range   {dem.min()} to {dem.max()}")
    if src.crs.is_geographic:
        lat = (src.bounds.bottom + src.bounds.top) / 2
        print(f"cell metres  {src.res[0] * 111_320 * math.cos(math.radians(lat)):.1f} x "
              f"{src.res[1] * 110_574:.1f}")
CRS          EPSG:4326 (geographic: True)
res          (0.00041666666666666664, 0.0002777777777777778)
nodata       None
dtype        int16
elev range   -32768 to 1074
cell metres  27.9 x 30.7

Three faults visible at once: degrees being used as a distance, an undeclared -32768 NoData, and an integer dtype.

Symptom Cause Fix
mean near 89°, max exactly 90° cell size in degrees convert to metres by latitude
ring of 90° around water or voids NoData in the 3×3 window mask before computing, dilate after
plausible but ~4% off square-cell assumption use both cell_x and cell_y
~15% too shallow via GDAL -s 111120 on a lat-adjusted DEM reproject first, or compute per-axis
slope fine, aspect rotated 90° np.gradient axis order rows are y, columns are x
huge phantom north-facing area flat cells arctan2(0, 0) is 0 — mask it
cliffs at every treeline it is a DSM, not a DTM use a bare-earth DTM
Seven slope symptoms mapped to their cause and fix, from degree cell sizes to a DSM being used instead of a DTM.
None of these is an error in the slope formula. All seven are input problems.

Step-by-step solution

1. Check the cell units before anything else 

if src.crs.is_geographic:
    raise ValueError(f"res {src.res} is in DEGREES — a gradient computed from it "
                     f"is ~90,000x too large")

This is the cause of the 88.79° result, and the arithmetic is worth seeing. A 100 m rise over one cell:

rise = 100.0
print(f"over 0.000417 degrees: gradient {rise / 0.000417:>12,.0f} -> "
      f"{math.degrees(math.atan(rise / 0.000417)):.4f}°")
print(f"over 27.9 metres:      gradient {rise / 27.9:>12,.2f} -> "
      f"{math.degrees(math.atan(rise / 27.9)):.4f}°")
over 0.000417 degrees: gradient      239,808 -> 89.9998°
over 27.9 metres:      gradient         3.58 -> 74.4272°

arctan of 239,808 is indistinguishable from 90°, so every cell with any relief at all saturates. That is why the mean is 88.79° rather than something merely too large.

2. Find the NoData that is not declared 

import numpy as np

print("declared:", src.nodata)
print("values below -100:", np.unique(dem[dem < -100])[:5])
print("count:", (dem < -100).sum())
declared: None
values below -100: [-32768]
count: 3,412

-32768 is a common int16 NoData sentinel. Left in place it produces a height difference of 33,000 m across a 28 m cell — a gradient of 1,200, which is 89.95°.

And it does not affect only those cells. Every 3×3 window containing one is contaminated, so a NoData lake gets a one-cell ring of 90° cliffs around its entire shoreline:

dem = np.where(dem < -100, np.nan, dem)

3. Mask the output one cell wider than the input 

from scipy.ndimage import binary_dilation

invalid = ~np.isfinite(dem)
halo = binary_dilation(invalid, np.ones((3, 3), bool))
slope = np.where(halo, np.nan, slope)

Masking only the original NoData cells leaves the ring behind, because those neighbouring cells had valid elevations — their windows did not. The dilation is what removes the halo.

4. Use both cell dimensions 

print(f"cell_x {cell_x:.1f} m, cell_y {cell_y:.1f} m, ratio {cell_x / cell_y:.3f}")
cell_x 27.9 m, cell_y 30.7 m, ratio 0.909

Copernicus DEM adjusts longitude spacing by latitude band so ground cells stay near-square — "near" being 9% out at 53°N. Using cell_x for both axes gives:

using both:      mean 21.33°  over 45°: 2.8%
using cell_x:    mean 22.28°  over 45°: 3.9%

A 4.5% error in the mean and a 39% error in "how much land is steeper than 45°". Nothing about the output looks wrong.

5. Check the aspect separately

Slope can be right while aspect is wrong, because slope uses the gradient magnitude and aspect uses its direction.

import pandas as pd

octants = pd.cut(aspect_valid, bins=np.arange(-22.5, 361, 45),
                 labels=["N", "NE", "E", "SE", "S", "SW", "W", "NW", "N2"], ordered=False)
print(octants.value_counts(normalize=True).sort_index().round(3).to_string())

A real landscape has an uneven but not extreme aspect distribution. Two red flags:

  • One octant far above 30% — usually flat cells defaulting to north.
  • The distribution rotated by 90° compared with what the terrain looks like — an axis-order mistake in np.gradient.
A NoData lake in a DEM producing a one-cell ring of ninety degree slopes around its shoreline.
The lake cells were masked. The ring around them was not — their windows touched the void.

Code examples

Example 1 — a diagnostic that names the fault 

import math

import numpy as np
import rasterio


def diagnose_slope(dem_path, slope=None, *, nodata_below=-100):
    problems = []

    with rasterio.open(dem_path) as src:
        raw = src.read(1)
        crs, res, declared, bounds = src.crs, src.res, src.nodata, src.bounds

    if crs is None:
        problems.append("no CRS — cell units unknown")
    elif crs.is_geographic:
        lat = (bounds.bottom + bounds.top) / 2
        cell_x = res[0] * 111_320 * math.cos(math.radians(lat))
        cell_y = res[1] * 110_574
        problems.append(
            f"geographic CRS: res {res} is DEGREES. Ground cell is "
            f"{cell_x:.1f} x {cell_y:.1f} m at {lat:.2f}°N"
        )
        if abs(cell_x / cell_y - 1) > 0.02:
            problems.append(f"cells are not square (ratio {cell_x / cell_y:.3f}) — "
                            f"use both dimensions, not one")
    else:
        cell_x, cell_y = res

    if np.issubdtype(raw.dtype, np.integer):
        problems.append(f"dtype is {raw.dtype} — cast to float64 before dividing")

    hidden = raw < nodata_below
    if declared is None and hidden.any():
        problems.append(
            f"nodata is None but {hidden.sum():,} cells are below {nodata_below} "
            f"(values {np.unique(raw[hidden])[:3]}) — undeclared NoData"
        )

    if slope is not None:
        finite = slope[np.isfinite(slope)]
        if finite.size:
            if finite.max() > 89.9:
                problems.append(f"max slope {finite.max():.2f}° — cell units are wrong")
            elif finite.max() > 80:
                problems.append(f"max slope {finite.max():.2f}° — implausible for real terrain")
            steep = (finite > 45).mean()
            if steep > 0.3:
                problems.append(f"{steep:.0%} of cells over 45° — implausible")
            print(f"  slope: mean {finite.mean():.2f}°  max {finite.max():.2f}°  "
                  f"over 45°: {steep:.1%}")

    print(f"  {crs} · res {res} · dtype {raw.dtype} · "
          f"elev {raw.min()} to {raw.max()}")
    for problem in problems:
        print(f"  ✗ {problem}")
    if not problems:
        print("  ✓ inputs look sound")
    return problems


diagnose_slope("dem.tif", slope)
  EPSG:4326 · res (0.0004166666666666666, 0.0002777777777777778) · dtype int16 · elev -32768 to 1074
  slope: mean 88.79°  max 90.00°  over 45%: 98.7%
  ✗ geographic CRS: res (0.0004166666666666666, 0.0002777777777777778) is DEGREES. Ground cell is 27.9 x 30.7 m at 53.07°N
  ✗ cells are not square (ratio 0.909) — use both dimensions, not one
  ✗ dtype is int16 — cast to float64 before dividing
  ✗ nodata is None but 3,412 cells are below -100 (values [-32768]) — undeclared NoData
  ✗ max slope 90.00° — cell units are wrong
  ✗ 99% of cells over 45° — implausible

Six findings from one call, four of them independent faults. Every one would have to be fixed before the slope raster meant anything.

Example 2 — removing the NoData halo 

from scipy.ndimage import binary_dilation


def masked_slope(dem, cell_x, cell_y, *, nodata_below=-100):
    invalid = ~np.isfinite(dem) | (dem < nodata_below)

    filled = np.where(invalid, np.nanmean(dem[~invalid]), dem)
    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:]
    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)
    slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))

    naive = np.where(invalid, np.nan, slope)
    halo = binary_dilation(invalid, np.ones((3, 3), bool))
    correct = np.where(halo, np.nan, slope)

    print(f"  masking only NoData:  max {np.nanmax(naive):.2f}°  "
          f"cells over 80°: {np.nansum(naive > 80):,}")
    print(f"  masking the halo too: max {np.nanmax(correct):.2f}°  "
          f"cells over 80°: {np.nansum(correct > 80):,}")
    print(f"  the halo cost {halo.sum() - invalid.sum():,} extra cells")
    return correct


slope = masked_slope(dem, cell_x=27.9, cell_y=30.7)
  masking only NoData:  max 84.42°  cells over 80°: 47
  masking the halo too: max 65.70°  cells over 80°: 0
  the halo cost 228 extra cells

Every cell over 80° disappears — all 47 of them were halo cells, and the maximum drops from 84.42° to 65.70°, which is the real terrain maximum. Losing 228 cells to NaN is a far better outcome than keeping 47 fake cliffs and a corrupted maximum.

The size of the effect depends on how the void boundary is oriented relative to the slope, which is why the maximum here is 84° rather than the 90° a deep -32768 sentinel would produce.

Example 3 — catching the aspect axis-order mistake 

def compare_aspect(dem, cell_x, cell_y):
    """np.gradient takes spacings in ARRAY order — rows (y) first."""
    dy_right, dx_right = np.gradient(dem, cell_y, cell_x)
    dy_wrong, dx_wrong = np.gradient(dem, cell_x, cell_y)      # swapped

    def summarise(dz_dx, dz_dy, label):
        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)
        counts, _ = np.histogram(aspect[np.isfinite(aspect)], bins=np.linspace(0, 360, 9))
        share = counts / counts.sum()
        names = ["N", "NE", "E", "SE", "S", "SW", "W", "NW"]
        print(f"  {label:16} slope mean {np.nanmean(slope):6.2f}°  aspect: "
              + " ".join(f"{n}{v:.0%}" for n, v in zip(names, share)))

    summarise(dx_right, dy_right, "correct order")
    summarise(dx_wrong, dy_wrong, "swapped")


compare_aspect(dem_clean, cell_x=27.9, cell_y=30.7)
  correct order    slope mean  21.55°  aspect: N8% NE11% E11% SE13% S13% SW15% W16% NW12%
  swapped          slope mean  21.76°  aspect: N7% NE12% E13% SE12% S12% SW16% W19% NW10%

The slopes differ by 0.21° — invisible. The aspect distribution shifts subtly and stays entirely plausible. This is why the axis-order mistake survives review: nothing about either output is obviously wrong, and the aspect map is quietly rotated.

The defence is not inspection, it is not using np.gradient for this. The explicit nine-slice form in the slope how-to has no axis ambiguity.

Explanation

Why the degree error saturates rather than merely exaggerating

If the cell size were merely 10× too small you would see slopes 10× too steep and notice immediately. The degree error is about 90,000× too small, and arctan compresses everything above a gradient of roughly 50 into the last degree before vertical.

So the output is not "too steep" in a way that scales — it is uniformly 89-point-something everywhere there is any relief at all, with a mean of 88.79° and a max of exactly 90.00°. That flatness at the top is the signature: a genuinely steep landscape has a spread of values, not a spike at the maximum.

Why the square-cell error is the dangerous one

Every other fault on the list is either obvious (90° everywhere) or localised (a ring around a lake). Assuming square cells produces a slope raster that is plausible everywhere and wrong by a few percent.

Where it bites is thresholds. A 4.5% shift in the mean becomes a 39% shift in "land steeper than 45°", because the threshold sits on the tail of the distribution where a small horizontal shift moves a lot of cells. Any analysis with a slope cutoff — buildable land, machine-accessible forestry, avalanche terrain — inherits that amplification.

Why a DSM produces cliffs at treelines

A DSM records the canopy top. At a forest edge the surface drops 15–25 m over one or two cells, which is a genuine gradient of the surface and a completely fictional gradient of the ground.

The signature is unmistakable once you look for it: slope maxima that trace field boundaries, hedgerows and plantation edges rather than following the terrain. If your steepest cells form long straight lines through gentle countryside, you have a DSM.

Global DEMs — Copernicus, SRTM, ASTER — are all DSMs. See DEM, DSM and DTM.

A 100 metre rise producing a gradient of 3.58 over 27.9 metres and 239,808 over 0.000417 degrees, saturating arctan at 90 degrees.
The signature is a spike at exactly 90°, not a distribution shifted upward.

Why to check against something you know

Every fault here is caught by one comparison against reality:

print(f"max elevation {np.nanmax(dem):.1f} m (Snowdon is 1085 m)")
print(f"max slope {np.nanmax(slope):.2f}° (real terrain rarely exceeds 70° over 30 m)")
print(f"land over 45°: {np.nanmean(slope > 45):.1%} (a few percent, even in mountains)")
max elevation 1074.5 m (Snowdon is 1085 m)
max slope 65.70° (real terrain rarely exceeds 70° over 30 m)
land over 45°: 2.8% (a few percent, even in mountains)

Three numbers, three known reference points. This catches unit errors, NoData contamination, datum problems and dtype truncation in one pass, and it takes longer to read than to run.

Edge cases or notes

  • A max of exactly 90.00° is always a bug. Real terrain over a 30 m cell tops out near 70°.
  • nodata=None does not mean there is none. Check for -32768, -9999 and 0-as-sea.
  • Integer DEMs truncate the gradient. Cast to float64 first.
  • The NoData halo is one cell wide for a 3×3 window. Dilate the mask before applying it.
  • np.gradient spacings are in array order — rows then columns. Swapping them rotates aspect and leaves slope essentially unchanged.
  • Flat cells have no aspect. Guard on the gradient magnitude, or you get a phantom north-facing plateau.
  • gdaldem -s cannot fix a latitude-adjusted geographic DEM — one factor for two different ground spacings. Reproject or compute per-axis.
  • Slope statistics are resolution-dependent. A 1 m DTM and a 30 m DEM give genuinely different answers for the same hillside; report the source resolution.

FAQ

Why is my maximum slope exactly 90°?

The cell size is in degrees. The gradient is roughly 90,000 times too large and arctan saturates. Convert the cell size to metres using the latitude.

Why is there a ring of cliffs around my lake?

NoData inside the 3×3 window. The lake cells were masked; their neighbours were not, but their windows touched the void. Dilate the invalid mask by one cell before applying it.

My slope looks fine but is a few percent off. What causes that?

Assuming square cells. On a Copernicus DEM at 53°N the ground cell is 27.9 × 30.7 m, and using one value for both shifts the mean by 4.5% and the "over 45°" share by 39%.

Why is my aspect map mostly north-facing?

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

The slope is right but the aspect looks rotated. Why?

np.gradient takes spacings in array order — rows (y) first. Swapping them rotates aspect while leaving slope almost unchanged, so it survives casual review.

Why do my steepest cells form straight lines across farmland?

You have a DSM, and those lines are hedgerows and field boundaries. Use a bare-earth DTM for terrain questions.

How do I know my slope raster is right?

Check three things against reality: the maximum elevation against a known summit, the maximum slope against about 70°, and the share of land over 45° against a few percent.