Hillshade Looks Flat, Washed Out or Inverted

Problem statement

The hillshade rendered, and it is not usable:

• Nearly uniform pale grey. There is relief in the DEM and none in the image.
• Hard black and white with no midtones, like a threshold rather than a shaded surface.
• Hills look like valleys. The relief is inverted, and once you see it you cannot un-see it.
• A ring of hard shadow around every lake and void.
• A texture of hedgerows and buildings where you expected landform.

Nothing errors. A hillshade always produces an image, and every one of these is an input or a parameter problem rather than a bug in the formula.

Quick answer

The standard deviation of the output diagnoses most of it in one number:

hs = hillshade(dem, 27.9, 30.7, altitude=80)          # the Snowdon massif, sun at 80°
finite = hs[np.isfinite(hs)]
print(f"mean {finite.mean():.1f}  sd {finite.std():.1f}  "
      f"range {finite.min():.0f}-{finite.max():.0f}")

mean 228.0  sd 25.9  range 98-255

A healthy hillshade has a standard deviation around 40–60 and spans most of 0–255. This one is 26 and never gets darker than 98 — washed out, because the sun altitude is too high.

Measure contrast, rule out saturation, check the azimuth

1. Measure the contrast rather than judging it by eye

import numpy as np


def contrast_report(hs, label=""):
    finite = hs[np.isfinite(hs)]
    dark = (finite < 20).mean()
    bright = (finite > 235).mean()
    verdict = ("washed out" if finite.std() < 30 else
               "saturated" if finite.std() > 75 else "healthy")
    print(f"  {label:20} mean {finite.mean():6.1f}  sd {finite.std():5.1f}  "
          f"range {finite.min():3.0f}-{finite.max():3.0f}  "
          f"dark {dark:5.2%}  bright {bright:5.2%}  -> {verdict}")


for altitude in (10, 30, 45, 60, 80):
    contrast_report(hillshade(dem, 27.9, 30.7, altitude=altitude), f"altitude {altitude}°")

  altitude 10°         mean   49.3  sd  49.5  range   0-231  dark 40.75%  bright 0.00%  -> healthy
  altitude 30°         mean  111.4  sd  61.2  range   0-254  dark 10.99%  bright 0.36%  -> healthy
  altitude 45°         mean  159.1  sd  55.3  range   0-255  dark 1.35%  bright 3.99%  -> healthy
  altitude 60°         mean  197.7  sd  43.1  range  13-255  dark 0.00%  bright 21.47%  -> healthy
  altitude 80°         mean  228.0  sd  25.9  range  98-255  dark 0.00%  bright 52.14%  -> washed out

Read the two clipping columns. At 80°, 52% of the image is above 235 — near-white, carrying no information. At 10°, 41% is below 20 — near-black, equally uninformative. Both extremes lose half the map; the difference is which end.

2. Rule out the saturated case

with rasterio.open("dem.tif") as src:
    res_x, res_y = src.res                                # degrees, not metres
p = np.pad(dem, 1, mode="edge")
dz_dx = ((p[:-2, 2:] + 2 * p[1:-1, 2:] + p[2:, 2:]) - (p[:-2, :-2] + 2 * p[1:-1, :-2] + p[2:, :-2])) / (8 * res_x)
dz_dy = ((p[2:, :-2] + 2 * p[2:, 1:-1] + p[2:, 2:]) - (p[:-2, :-2] + 2 * p[:-2, 1:-1] + p[:-2, 2:])) / (8 * res_y)
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
print(f"slope max {slope.max():.2f}°, over 45°: {(slope > 45).mean():.1%}")
contrast_report(hillshade(dem, res_x, res_y), "cell size in degrees")

slope max 90.00°, over 45°: 98.7%
  cell size in degrees mean   58.0  sd  70.5  range   0-182  dark 54.30%  bright 0.00%  -> healthy

If the slopes have saturated at 90°, the hillshade becomes nearly binary — a cell is either lit or in full shadow, with little between. Note what the report says: the standard deviation is not extreme, and the verdict is "healthy". The tell is the dark column — half the map is pure black on a mountain lit at 45°.

The cause is always the same: the cell size is in degrees. See slope values are wrong.

3. Check the azimuth for relief inversion

This one cannot be diagnosed from a number — it is a perceptual failure and the statistics are perfectly healthy.

for azimuth in (315, 135):
    contrast_report(hillshade(dem, 27.9, 30.7, azimuth=azimuth), f"azimuth {azimuth}°")

  azimuth 315°         mean  159.1  sd  55.3  range   0-255  dark 1.35%  bright 3.99%  -> healthy
  azimuth 135°         mean  170.4  sd  50.2  range   0-255  dark 0.54%  bright 10.36%  -> healthy

Both healthy. The 135° version is a technically correct rendering in which every hill reads as a hollow, because human vision assumes light from above-left.

There is no numeric test. The defence is procedural: use 315°, and if a particular ridge is badly lit, change the altitude or blend several directions rather than moving the light south.

4. Distinguish "flat image" from "flat terrain"

The Fens east of Chatteris, from the Environment Agency's lidar terrain model at 25 m:

with rasterio.open("fens_dtm.tif") as src:
    flat_dem = src.read(1, masked=True).filled(np.nan).astype("float64")
fen_slope = np.degrees(np.arctan(np.hypot(*np.gradient(flat_dem, 25, 25))))
print(f"elevation range {np.nanmax(flat_dem) - np.nanmin(flat_dem):.1f} m "
      f"over {flat_dem.shape[1] * 25 / 1000:.0f} km")
print(f"slope: mean {np.nanmean(fen_slope):.2f}°  p95 {np.nanpercentile(fen_slope, 95):.2f}°")

elevation range 36.0 m over 10 km
slope: mean 0.56°  p95 2.21°

If the 95th-percentile slope is under about 2°, the terrain really is flat and no amount of parameter tuning will produce relief. That is what z_factor is for:

for z in (1, 5, 10, 20):
    contrast_report(hillshade(flat_dem, 25, 25, z_factor=z), f"z_factor {z}")

  z_factor 1           mean  180.4  sd   2.3  range 125-207  dark 0.00%  bright 0.00%  -> washed out
  z_factor 5           mean  180.3  sd  11.0  range   0-251  dark 0.01%  bright 0.06%  -> washed out
  z_factor 10          mean  179.5  sd  20.3  range   0-255  dark 0.14%  bright 0.80%  -> washed out
  z_factor 20          mean  177.0  sd  32.9  range   0-255  dark 0.95%  bright 2.83%  -> healthy

Twenty-fold exaggeration on a floodplain is legitimate cartography. It is not a measurement — never compute slope from a z-exaggerated DEM.

5. Look for the NoData ring

Copernicus has no voids, so take a copy with the lakes under Snowdon blanked out, as many DEMs deliver them:

from scipy.ndimage import binary_dilation

lakes = np.where(water == 80, np.nan, dem)
hs = hillshade(np.where(np.isnan(lakes), np.nanmean(lakes), lakes), 27.9, 30.7)
finite = hs[np.isfinite(lakes)]
invalid = ~np.isfinite(lakes)
halo = binary_dilation(invalid, np.ones((3, 3), bool)) & ~invalid
print(f"{halo.sum():,} halo cells; their hillshade sd is "
      f"{hs[halo][np.isfinite(hs[halo])].std():.1f} vs {finite.std():.1f} overall")

600 halo cells; their hillshade sd is 86.1 vs 56.3 overall

The ring's statistics differ from the rest of the map because each of its cells is shaded against a fake surface — here the mean elevation poured into the lakes, elsewhere a -32768 — rather than the shoreline. Dilate the invalid mask by one cell before applying it.

Snowdon and the Fens: a diagnostic and an altitude search

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, read from its public bucket; ESA WorldCover's water class on the same grid, to blank the lakes for step 5; and the Environment Agency's lidar terrain model of the Fens east of Chatteris at 25 m, for genuinely flat ground. hillshade() is the Horn-based function from the hillshade how-to, with a z_factor.

import math
import numpy as np
import rasterio
import requests
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)
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)
dem = mosaic[0].astype("float64")

water = np.zeros(dem.shape, 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)

r = requests.get(
    "https://environment.data.gov.uk/spatialdata/lidar-composite-digital-terrain-model-dtm-1m/wcs",
    params=[("service", "WCS"), ("version", "2.0.1"), ("request", "GetCoverage"),
            ("CoverageId", "13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m"),
            ("format", "image/tiff"), ("subset", "E(540000,550000)"), ("subset", "N(280000,290000)"),
            ("SCALEFACTOR", "0.04")], timeout=300)
r.raise_for_status()
with open("fens_dtm.tif", "wb") as f:
    f.write(r.content)


def hillshade(dem, cell_x, cell_y, *, azimuth=315, altitude=45, z_factor=1.0):
    p = np.pad(dem * z_factor, 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.arctan(np.hypot(dz_dx, dz_dy))
    aspect = np.arctan2(dz_dy, -dz_dx)
    az, alt = math.radians(360.0 - azimuth + 90.0), math.radians(altitude)
    shaded = math.sin(alt) * np.cos(slope) + math.cos(alt) * np.sin(slope) * np.cos(az - aspect)
    return np.clip(shaded, 0, 1) * 255


print(f"dem.tif {dem.shape[1]} × {dem.shape[0]}, {dem.min():.0f}–{dem.max():.0f} m; "
      f"{(water == 80).sum():,} lake cells; fens_dtm.tif downloaded")

dem.tif 240 × 252, 54–1075 m; 1,327 lake cells; fens_dtm.tif downloaded

Example 1 — a diagnostic that names the fault

import math

import numpy as np
import rasterio


def diagnose_hillshade(dem_path, hs, *, azimuth=315, altitude=45, z_factor=1.0):
    problems = []

    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype("float64")
        if src.nodata is not None:
            dem = np.where(dem == src.nodata, np.nan, dem)
        geographic = src.crs.is_geographic
        res = src.res
        lat = (src.bounds.bottom + src.bounds.top) / 2

    if geographic:
        cell_x = res[0] * 111_320 * math.cos(math.radians(lat))
        cell_y = res[1] * 110_574
        problems.append(f"geographic CRS — cell must be {cell_x:.1f} x {cell_y:.1f} m, "
                        f"not {res}")

    finite = hs[np.isfinite(hs)]
    sd, mean = finite.std(), finite.mean()
    dark, bright = (finite < 20).mean(), (finite > 235).mean()

    black = (finite == 0).mean()
    if black > 0.3:
        problems.append(f"{black:.0%} of cells pure black — slopes have saturated; check the cell units")
    elif sd < 30:
        relief = np.nanmax(dem) - np.nanmin(dem)
        if relief < 50:
            problems.append(f"sd {sd:.1f} with only {relief:.0f} m of relief — "
                            f"terrain is genuinely flat; raise z_factor")
        elif altitude > 55:
            problems.append(f"sd {sd:.1f} with altitude {altitude}° — sun too high; "
                            f"drop to 30-45°")
        else:
            problems.append(f"sd {sd:.1f} — low contrast; try z_factor or a lower altitude")

    if bright > 0.4:
        problems.append(f"{bright:.0%} of cells above 235 — washed out")
    if dark > 0.4 and black <= 0.3:
        problems.append(f"{dark:.0%} of cells below 20 — mostly in shadow; raise the altitude")

    if not (270 <= azimuth <= 360 or azimuth == 0):
        problems.append(f"azimuth {azimuth}° is outside the 270-360° band — "
                        f"expect relief inversion; the convention is 315°")

    if z_factor != 1.0:
        problems.append(f"z_factor {z_factor} — display only, never measure from this")

    print(f"  az {azimuth}° alt {altitude}° · mean {mean:.1f} sd {sd:.1f} "
          f"dark {dark:.1%} bright {bright:.1%}")
    for problem in problems:
        print(f"  ✗ {problem}")
    if not problems:
        print("  ✓ hillshade parameters look sound")
    return problems


diagnose_hillshade("dem.tif", hillshade(dem, 27.9, 30.7, altitude=80), altitude=80)

  az 315° alt 80° · mean 228.0 sd 25.9 dark 0.0% bright 52.1%
  ✗ geographic CRS — cell must be 27.9 x 30.7 m, not (0.00041666666666666664, 0.0002777777777777778)
  ✗ sd 25.9 with altitude 80° — sun too high; drop to 30-45°
  ✗ 52% of cells above 235 — washed out

The altitude branch is what turns "the image looks pale" into a specific instruction. The first finding is a standing reminder rather than this image's fault — the hillshade was computed with the metre cell sizes — but a diagnostic that only sees the file cannot know that, so it says so.

Example 2 — finding a good altitude automatically

def best_altitude(dem, cell_x, cell_y, *, azimuth=315, candidates=(20, 25, 30, 35, 40, 45, 50, 55)):
    """Pick the altitude that maximises contrast without losing too much to shadow."""
    import pandas as pd

    rows = []
    for altitude in candidates:
        hs = hillshade(dem, cell_x, cell_y, azimuth=azimuth, altitude=altitude)
        finite = hs[np.isfinite(hs)]
        dark, bright = (finite < 20).mean(), (finite > 235).mean()
        # penalise clipping at both ends — those cells carry no information
        score = finite.std() * (1 - dark - bright)
        rows.append({"altitude": altitude, "sd": round(finite.std(), 1),
                     "dark": f"{dark:.1%}", "bright": f"{bright:.1%}",
                     "score": round(score, 1)})

    frame = pd.DataFrame(rows)
    best = frame.loc[frame["score"].idxmax(), "altitude"]
    print(frame.to_string(index=False))
    print(f"\nbest altitude: {best}°")
    return int(best)


best_altitude(dem, 27.9, 30.7)

 altitude   sd  dark bright  score
       20 58.2 23.3%   0.1%   44.6
       25 60.4 16.5%   0.2%   50.4
       30 61.2 11.0%   0.4%   54.3
       35 60.4  6.4%   0.8%   56.1
       40 58.4  3.4%   1.8%   55.4
       45 55.3  1.4%   4.0%   52.4
       50 51.6  0.4%   7.9%   47.3
       55 47.4  0.1%  13.6%   40.9

best altitude: 35°

Thirty-five degrees for this terrain, against the conventional 45°. Note that raw standard deviation peaks at 30° while the score peaks at 35°: at 30° an extra 4.6% of the map has clipped to black, and those cells contribute to the spread without carrying any shape. The score penalises clipped cells at both ends, because a cell at 0 or 255 carries no shape information regardless of how much it contributes to the standard deviation.

Use this as a starting point rather than a rule — for a map that will sit alongside others, consistency with 45° may matter more than optimal contrast.

Example 3 — proving the inversion is real

def inversion_demo(dem, cell_x, cell_y, ridge_row=120):
    """A single ridge profile, shaded from both directions."""
    for azimuth, label in [(315, "NW — reads as a ridge"), (135, "SE — reads as a valley")]:
        hs = hillshade(dem, cell_x, cell_y, azimuth=azimuth)
        profile = hs[ridge_row, 100:140]
        elev = dem[ridge_row, 100:140]
        peak = int(np.argmax(elev))
        print(f"  {label}")
        print(f"    elevation peaks at column {peak} ({elev[peak]:.0f} m)")
        print(f"    hillshade there is {profile[peak]:.0f}; "
              f"west side {profile[max(0, peak - 6)]:.0f}, "
              f"east side {profile[min(len(profile) - 1, peak + 6)]:.0f}")


inversion_demo(dem, 27.9, 30.7)

  NW — reads as a ridge
    elevation peaks at column 36 (504 m)
    hillshade there is 206; west side 223, east side 157
  SE — reads as a valley
    elevation peaks at column 36 (504 m)
    hillshade there is 146; west side 123, east side 199

The flanks are mirrored. Under north-west lighting the west face is bright (223) and the east face dark (157), which the eye reads as a ridge. Under south-east lighting the pattern reverses — west 123, east 199 — and the identical landform reads as a hollow.

The elevation is unchanged. Only which flank is lit swaps, and that is enough to invert the perception entirely.

Why sun altitude controls contrast

Why sun altitude controls contrast

The hillshade formula has two terms:

sin(altitude) * cos(slope)                            # ambient, same everywhere
+ cos(altitude) * sin(slope) * cos(azimuth - aspect)  # directional

At 80° altitude, sin(80°) is 0.98 and cos(80°) is 0.17 — the ambient term dominates and the directional term, which carries all the shape information, is scaled down by a factor of six. Everything is lit and nothing is distinguished.

At 10°, cos(10°) is 0.98 and the directional term dominates, but sin(10°) is 0.17, so the base level is very dark and any slope facing away from the sun clips to zero.

Between roughly 30° and 45° both terms are substantial. The measured optimum on this terrain was 35°, with the conventional 45° close behind.

Why relief inversion has no numeric signature

Both azimuths produce statistically identical images — the same standard deviation, the same range, the same distribution shape. The difference is entirely in which flank is bright.

Human vision resolves shape from shading using a hard-wired assumption that illumination comes from above. In an image, "above" is the top of the frame, and the assumption extends slightly leftward. North-west lighting matches it; south-east lighting contradicts it, and the visual system resolves the contradiction by flipping the interpreted surface.

This is why 315° is a requirement rather than a preference, and why no amount of measurement will detect the failure.

Why a hillshade of a DSM looks like fabric

A DSM records the top of whatever is there. Over farmland that means every hedgerow is a 3–5 m wall and every building a block. In a hillshade those become sharp bright-dark pairs, thousands of them, at a scale much finer than the landform.

The visual result is a woven texture that dominates the image. The landform is still present underneath, and it is unreadable.

The signature: bright-dark pairs following straight lines and field boundaries rather than contours. If your hillshade looks like corduroy, you have a DSM — see DEM, DSM and DTM.

Why hillshade is the best DEM quality check

Shading is a derivative operation, so it amplifies anything with a sharp edge. That makes it useless for measurement and superb for inspection:

• Tile seams from mosaicking appear as straight lines.
• Interpolated voids appear as unnaturally smooth patches.
• Striping from sensor artefacts becomes obvious.
• Resampling artefacts appear as a faint grid.

None of these is visible in an elevation colour ramp, where a 2 m step in a 1,000 m range is one part in five hundred. In a hillshade the same step is a hard bright line. Render a hillshade of every DEM before trusting it.

Healthy spread, z_factor and two-spike histograms

• Standard deviation 40–60 is healthy and under 30 is washed out — but saturated slopes can hide inside a normal-looking standard deviation (70.5 here). Check the share of pure-black cells too.
• Relief inversion has no numeric test. Keep the azimuth in the 270–360° band.
• z_factor is for display only. Never derive slope, aspect or volume from an exaggerated DEM.
• The NoData halo has a much higher local standard deviation than the image overall — a useful confirmation.
• clip(0, 1) before scaling stops cells facing directly away from the sun wrapping to bright.
• A histogram with two spikes at 0 and 255 means saturation, not contrast.
• Optimal altitude is terrain-dependent — 30° here, higher for gentler ground. Consistency across a map series may matter more.
• Multidirectional blending reduces overall contrast but rescues detail on shadowed faces.

Internal links

• How to make a hillshade from a DEM in Python — the implementation with these guards
• Slope values are wrong or absurdly steep — the saturated-slope cause
• How to calculate slope and aspect from a DEM in Python — the two inputs
• Digital elevation models explained — why a DSM hillshade looks like fabric
• Slope and aspect explained — the aspect convention behind the azimuth
• Rasterio output is black or empty — a related "the image is wrong" problem
• The raster data model explained — NoData and dtype in the output
• Choropleth classification explained — the same contrast question in colour

FAQ

Why is my hillshade pale and featureless?

The sun altitude is too high. At 80° the standard deviation drops to 26 and 52% of cells are near-white. Use 30–45°.

Why is it hard black and white with no midtones?

The slopes have saturated at 90°, which means the cell size is in degrees rather than metres. Fix the cell size and the midtones return.

Why do the hills look like valleys?

The azimuth is in the southern half. Human vision assumes light from above-left, so south-east lighting inverts the perceived relief. Use 315°.

The statistics look fine but the image looks wrong. What now?

Check the azimuth. Relief inversion produces statistically identical output — there is no numeric test for it.

My terrain is genuinely flat. Can I still get a hillshade?

Yes, with z_factor. A floodplain may need 10–20× exaggeration. It is a visualisation, not a measurement.

Why is there a hard ring around my lakes?

NoData inside the 3×3 window. Those halo cells have a much higher local standard deviation than the rest of the image. Dilate the invalid mask by one cell.

What altitude is best?

It depends on the terrain — 35° scored best on the example here, against the conventional 45°. Score candidates by standard deviation penalised for clipped cells, and weigh consistency with other maps.