Slope Values Are Wrong or Absurdly Steep

Problem statement

The slope raster came out and the numbers cannot be right. This is the Snowdon massif, in a DEM exported the way many are — int16, lakes written as -32768, no NoData tag:

import numpy as np
import rasterio

with rasterio.open("dem.tif") as src:
    dy, dx = np.gradient(src.read(1).astype("float64"), *src.res[::-1])
naive_slope = np.degrees(np.arctan(np.hypot(dx, dy)))
print(f"mean {naive_slope.mean():.2f}°  max {naive_slope.max():.2f}°  over 45°: {(naive_slope > 45).mean():.1%}")

mean 88.34°  max 90.00°  over 45°: 98.2%

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 1075
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.

Cell units, hidden NoData, wider masks, both dimensions

1. Check the cell units before anything else

if src.crs.is_geographic:
    raise ValueError(f"res {src.res} is in DEGREES — convert to ground metres first")

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.4109°

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

print("declared:", src.nodata)
print("values below -100:", np.unique(dem[dem < -100])[:5])
print(f"count: {(dem < -100).sum():,}")

declared: None
values below -100: [-32768]
count: 1,327

-32768 is the usual int16 NoData sentinel — SRTM uses it for voids. Here it marks the lakes: Llyn Llydaw, Glaslyn and the others under Snowdon. Left in place, a lake cell sits 33,000 m below its shore across a 28 m cell — a gradient of about 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_clean = np.where(dem < -100, np.nan, dem.astype("float64"))

3. Mask the output one cell wider than the input

from scipy.ndimage import binary_dilation


def horn_slope(z, cell_x, cell_y):
    p = np.pad(z, 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 np.degrees(np.arctan(np.hypot(dz_dx, dz_dy))), dz_dx, dz_dy


invalid = ~np.isfinite(dem_clean)
filled = np.where(invalid, np.nanmean(dem_clean), dem_clean)     # keep the window arithmetic finite
slope, dz_dx, dz_dy = horn_slope(filled, 27.9, 30.7)
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

import math

lat = (src.bounds.bottom + src.bounds.top) / 2
cell_x = src.res[0] * 111_320 * math.cos(math.radians(lat))
cell_y = src.res[1] * 110_574
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.907

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:

for label, cy in [("using both", cell_y), ("using cell_x", cell_x)]:
    s = np.where(halo, np.nan, horn_slope(filled, cell_x, cy)[0])
    print(f"{label:14} mean {np.nanmean(s):.2f}°  over 45°: {np.nanmean(s > 45):.1%}")

using both     mean 21.84°  over 45°: 2.8%
using cell_x   mean 22.82°  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

gradient = np.hypot(dz_dx, dz_dy)
aspect = np.degrees(np.arctan2(dz_dy, -dz_dx)) % 360
aspect_valid = aspect[~halo & (gradient > 1e-8)]
octants = pd.Series(pd.cut(aspect_valid, bins=np.arange(-22.5, 383, 45),
                           labels=["N", "NE", "E", "SE", "S", "SW", "W", "NW", "N"], ordered=False))
print(octants.value_counts(normalize=True).reindex(["N", "NE", "E", "SE", "S", "SW", "W", "NW"]).round(3).to_string())

N     0.091
NE    0.099
E     0.104
SE    0.125
S     0.136
SW    0.135
W     0.142
NW    0.168

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 diagnostic for an int16 DEM with undeclared voids

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, read from its public bucket, then exported the way a lot of DEMs arrive: as int16, with every cell ESA WorldCover maps as open water set to -32768, and no NoData tag. The faults on this page are that export's; the terrain is real.

import numpy as np
import rasterio
from rasterio.merge import merge
from rasterio.warp import reproject, Resampling

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)
water = np.zeros(mosaic.shape[1:], dtype="uint8")
with rasterio.open("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
                   "ESA_WorldCover_10m_2021_v200_N51W006_Map.tif") as src:
    reproject(rasterio.band(src, 1), water, dst_transform=transform, dst_crs="EPSG:4326",
              resampling=Resampling.mode)
exported = np.where(water == 80, -32768, np.round(mosaic[0])).astype("int16")
with rasterio.open(tiles[0]) as src:
    profile = src.profile | {"driver": "GTiff", "dtype": "int16", "nodata": None, "count": 1,
                             "width": exported.shape[1], "height": exported.shape[0],
                             "transform": transform, "compress": "deflate"}
with rasterio.open("dem.tif", "w", **profile) as dst:
    dst.write(exported, 1)
print(f"dem.tif {exported.shape[1]} × {exported.shape[0]} int16, "
      f"{(exported == -32768).sum():,} lake cells written as -32768")

dem.tif 240 × 252 int16, 1,327 lake cells written as -32768

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", naive_slope)

  slope: mean 88.34°  max 90.00°  over 45°: 98.2%
  EPSG:4326 · res (0.00041666666666666664, 0.0002777777777777778) · dtype int16 · elev -32768 to 1075
  ✗ geographic CRS: res (0.00041666666666666664, 0.0002777777777777778) is DEGREES. Ground cell is 27.9 x 30.7 m at 53.06°N
  ✗ cells are not square (ratio 0.907) — use both dimensions, not one
  ✗ dtype is int16 — cast to float64 before dividing
  ✗ nodata is None but 1,327 cells are below -100 (values [-32768]) — undeclared NoData
  ✗ max slope 90.00° — cell units are wrong
  ✗ 98% 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.astype("float64"), cell_x=27.9, cell_y=30.7)

  masking only NoData:  max 82.21°  cells over 80°: 61
  masking the halo too: max 65.73°  cells over 80°: 0
  the halo cost 600 extra cells

Every cell over 80° disappears — all 61 of them were halo cells around the lakes, and the maximum drops from 82.21° to 65.73°, the real terrain maximum. Losing 600 cells to NaN is a far better outcome than keeping 61 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 82° rather than the 90° a raw -32768 would produce: the function filled the lakes with the mean elevation before the window arithmetic, so the halo is a drop to about 500 m, not to −32,768.

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  22.04°  aspect: N9% NE11% E11% SE13% S14% SW14% W17% NW12%
  swapped          slope mean  22.25°  aspect: N8% NE12% E13% SE11% 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.

Why degree cell sizes saturate slope

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 67,000 to 110,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.34° 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 by the height of the trees over one or two cells — the Copernicus DSM sits a median 10 m above the ground under Epping Forest's trees — 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.

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 1075.0 m (Snowdon is 1085 m)
max slope 65.73° (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.

A 90° maximum, NoData halos and DSM treelines

• 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.

Internal links

• How to calculate slope and aspect from a DEM in Python — the implementation with these guards in place
• Slope and aspect explained — what the four decisions are
• Digital elevation models explained — DSM cliffs at treelines
• Vertical datums explained — why slope is unaffected by the vertical zero
• The raster data model explained — dtype, NoData and the transform
• Hillshade looks flat, washed out or inverted — the same inputs, a different symptom
• Rasterio output is black or empty — NoData contaminating a different product
• Projected vs geographic CRS explained — why degrees are not distances

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.