⚙️ How-To Terrain & elevation Raster & rasterio Terrain & Elevation · step 12 of 12 Updated

How to Calculate a Viewshed in Python

Problem statement

A viewshed answers "what can be seen from here" — the terrain visible from an observer, given the shape of the ground between. It is the basis of turbine visibility assessments, telecoms coverage, and any planning question involving sightlines.

There is no rasterio.viewshed(). And the calculation has parameters that change the answer by a factor of two:

observer at the summit, 1,074.5 m
  eye height 0 m                visible 21.44%
  eye height 1.7 m              visible 23.72%
  30 m mast                     visible 38.64%

A 30 m mast sees 63% more of the landscape than a person standing on the same spot. That is a real physical effect and it is entirely invisible in the output raster — two viewsheds of the same point, produced by the same code, differing only in a parameter nobody recorded.

Quick answer

Cast a ray to every boundary cell and track the maximum vertical angle seen so far along it:

import numpy as np
from scipy.ndimage import map_coordinates

EARTH_RADIUS_M = 6_371_000.0


def viewshed(dem, row, col, cell_x, cell_y, *, observer_h=1.7, target_h=0.0,
             max_range=None, curvature=True, refraction=0.13):
    rows, cols = dem.shape
    visible = np.zeros(dem.shape, bool)
    visible[row, col] = True
    eye = dem[row, col] + observer_h

    boundary = ([(0, c) for c in range(cols)] + [(rows - 1, c) for c in range(cols)]
                + [(r, 0) for r in range(rows)] + [(r, cols - 1) for r in range(rows)])

    for br, bc in boundary:
        n = int(max(abs(br - row), abs(bc - col)))
        if n == 0:
            continue
        rr = np.linspace(row, br, n + 1)
        cc = np.linspace(col, bc, n + 1)
        d = np.hypot((rr - row) * cell_y, (cc - col) * cell_x)
        z = map_coordinates(dem, [rr, cc], order=1, mode="nearest")

        drop = (1 - refraction) * d ** 2 / (2 * EARTH_RADIUS_M) if curvature else 0.0
        angle = np.divide(z + target_h - drop - eye, d,
                          out=np.full_like(d, -np.inf), where=d > 0)
        if max_range:
            angle = np.where(d <= max_range, angle, -np.inf)

        horizon = np.maximum.accumulate(np.concatenate([[-np.inf], angle[:-1]]))
        seen = angle >= horizon

        ri = np.rint(rr).astype(int).clip(0, rows - 1)
        ci = np.rint(cc).astype(int).clip(0, cols - 1)
        visible[ri[seen], ci[seen]] = True

    return visible


vis = viewshed(dem, row, col, 27.9, 30.7)
print(f"visible {vis.mean():.2%} ({vis.sum():,} cells)")
visible 23.72% (14,285 cells)
Parameter Typical Effect
observer height 1.7 m (a person) large — 21.4% to 23.7% for the first 1.7 m
target height 0 m (the ground) large for tall targets
max range none, or a limit 5 km limit: 23.7% → 17.2%
curvature on negligible under ~10 km, essential beyond
A ray from an observer across terrain, with the running maximum angle forming a horizon that hides cells behind ridges.
A cell is visible if its angle above the observer exceeds every angle seen closer along the same ray.

Step-by-step solution

1. Understand the angle test

Visibility is not about height, it is about angle. A cell is visible if the vertical angle from the observer's eye to that cell exceeds the largest angle already encountered along the ray:

angle = (z + target_h - eye) / distance
horizon = np.maximum.accumulate(angle_so_far)
seen = angle >= horizon

That is why a distant mountain top is visible over a near hill, and a valley floor just beyond the same hill is not. The running maximum is the horizon as it appears from the observer, and it only ever rises along a ray.

np.maximum.accumulate computes it in one vectorised pass, which is what makes a pure-NumPy viewshed fast enough to be practical.

2. Set the observer height deliberately 

for h, label in [(0, "ground level"), (1.7, "a person"), (10, "a house"), (30, "a mast")]:
    vis = viewshed(dem, row, col, 27.9, 30.7, observer_h=h)
    print(f"  {label:14} {h:4} m: visible {vis.mean():6.2%} ({vis.sum():,} cells)")
  ground level      0 m: visible 21.44% (12,911 cells)
  a person        1.7 m: visible 23.72% (14,285 cells)
  a house          10 m: visible 33.34% (20,078 cells)
  a mast           30 m: visible 38.64% (23,270 cells)

The first 1.7 m adds 1,374 cells — 11% more visible area — because it lifts the eye above the local micro-relief that would otherwise block low-angle rays.

A default of 0 is almost never what anyone means. 1.7 m is the standard for a standing person; use the actual structure height for infrastructure.

3. Set the target height for what you are looking at

target_h raises every cell before testing it, which answers a different question: not "can I see the ground there" but "can I see a 100 m turbine standing there".

turbine_visibility = viewshed(dem, row, col, 27.9, 30.7,
                              observer_h=1.7, target_h=100.0)

For a wind farm assessment this is the parameter that matters, and the two viewsheds look completely different. Confusing them produces an assessment answering the wrong question.

4. Decide whether curvature matters 

for curvature in (True, False):
    vis = viewshed(dem, row, col, 27.9, 30.7, curvature=curvature)
    print(f"  curvature {str(curvature):5}: visible {vis.mean():6.2%}")
  curvature True : visible 23.72%
  curvature False: visible 23.77%

A difference of 0.05 percentage points — negligible, because this DEM spans about 8 km. The earth drops roughly d² / (2R) below a tangent plane: 5 m at 8 km, but 78 m at 32 km and 314 m at 64 km.

So: under about 10 km it does not matter; beyond 20 km it dominates. Refraction bends light back down and recovers roughly 13% of the drop, which is the (1 - 0.13) factor.

5. Limit the range if the question has one 

for limit in (None, 20_000, 10_000, 5_000):
    vis = viewshed(dem, row, col, 27.9, 30.7, max_range=limit)
    label = "unlimited" if limit is None else f"{limit / 1000:.0f} km"
    print(f"  {label:>9}: visible {vis.mean():6.2%}")
  unlimited: visible 23.72%
      20 km: visible 23.72%
      10 km: visible 23.72%
       5 km: visible 17.24%

The 5 km limit removes a quarter of the visible area. The 10 and 20 km limits change nothing, because the DEM extent is smaller than either.

A range limit is a modelling decision about what "visible" means — a turbine 30 km away is geometrically visible and may be perceptually irrelevant. State it.

Visible area rising from 21.4 percent at ground level to 38.6 percent from a 30 metre mast at the same location.
Same point, same terrain. Observer height alone changes the answer by 80% between ground level and a mast.

Code examples

Example 1 — a viewshed with every parameter recorded 

import math

import numpy as np
import rasterio
from rasterio.transform import rowcol
from scipy.ndimage import map_coordinates

EARTH_RADIUS_M = 6_371_000.0


def viewshed_from_dem(path, observer_xy, *, observer_h=1.7, target_h=0.0,
                      max_range=None, curvature=True, refraction=0.13):
    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)
        transform, crs = src.transform, src.crs
        if crs.is_geographic:
            lat = (src.bounds.bottom + src.bounds.top) / 2
            cell_x = abs(transform.a) * 111_320 * math.cos(math.radians(lat))
            cell_y = abs(transform.e) * 110_574
        else:
            cell_x, cell_y = abs(transform.a), abs(transform.e)
        profile = src.profile.copy()

    if np.isnan(dem).any():
        raise ValueError("voids block rays — fill the DEM before computing a viewshed")

    row, col = rowcol(transform, observer_xy[0], observer_xy[1])
    rows, cols = dem.shape
    if not (0 <= row < rows and 0 <= col < cols):
        raise ValueError(f"observer at {observer_xy} is outside the DEM")

    visible = np.zeros(dem.shape, bool)
    visible[row, col] = True
    eye = dem[row, col] + observer_h

    boundary = ([(0, c) for c in range(cols)] + [(rows - 1, c) for c in range(cols)]
                + [(r, 0) for r in range(rows)] + [(r, cols - 1) for r in range(rows)])

    for br, bc in boundary:
        n = int(max(abs(br - row), abs(bc - col)))
        if n == 0:
            continue
        rr, cc = np.linspace(row, br, n + 1), np.linspace(col, bc, n + 1)
        d = np.hypot((rr - row) * cell_y, (cc - col) * cell_x)
        z = map_coordinates(dem, [rr, cc], order=1, mode="nearest")

        drop = (1 - refraction) * d ** 2 / (2 * EARTH_RADIUS_M) if curvature else 0.0
        angle = np.divide(z + target_h - drop - eye, d,
                          out=np.full_like(d, -np.inf), where=d > 0)
        if max_range:
            angle = np.where(d <= max_range, angle, -np.inf)

        horizon = np.maximum.accumulate(np.concatenate([[-np.inf], angle[:-1]]))
        seen = angle >= horizon
        ri = np.rint(rr).astype(int).clip(0, rows - 1)
        ci = np.rint(cc).astype(int).clip(0, cols - 1)
        visible[ri[seen], ci[seen]] = True

    cell_km2 = cell_x * cell_y / 1e6
    meta = {
        "observer_x": observer_xy[0], "observer_y": observer_xy[1],
        "observer_elev_m": round(float(dem[row, col]), 2),
        "observer_h_m": observer_h, "target_h_m": target_h,
        "max_range_m": max_range or "unlimited",
        "curvature": curvature, "refraction": refraction,
        "visible_cells": int(visible.sum()),
        "visible_km2": round(float(visible.sum() * cell_km2), 2),
        "visible_share": round(float(visible.mean()), 4),
    }
    print(f"  observer {dem[row, col]:.1f} m + {observer_h} m, target +{target_h} m")
    print(f"  visible {meta['visible_share']:.2%} = {meta['visible_km2']} km² "
          f"({meta['visible_cells']:,} cells)")
    return visible, meta, profile


summit = (-4.0763, 53.0685)
vis, meta, profile = viewshed_from_dem("snowdonia_glo30.tif", summit)
  observer 1074.5 m + 1.7 m, target +0.0 m
  visible 23.72% = 12.24 km² (14,285 cells)

The DEM covers 51.6 km², so just under a quarter of it is visible from the summit.

The metadata dictionary is the deliverable as much as the raster is. A viewshed without its observer height, target height and range limit cannot be compared with any other viewshed or defended in a planning objection.

Example 2 — how sensitive the answer is 

import pandas as pd


def viewshed_sensitivity(path, observer_xy):
    rows = []
    for observer_h in (0, 1.7, 10, 30):
        for max_range in (None, 5_000):
            _, meta, _ = viewshed_from_dem(path, observer_xy,
                                           observer_h=observer_h, max_range=max_range)
            rows.append({
                "observer_h": observer_h,
                "range": "unlimited" if max_range is None else f"{max_range // 1000} km",
                "visible": f"{meta['visible_share']:.2%}",
                "km2": meta["visible_km2"],
            })
    frame = pd.DataFrame(rows).pivot(index="observer_h", columns="range", values="km2")
    print(frame.to_string())
    return frame


viewshed_sensitivity("snowdonia_glo30.tif", summit)
range        5 km  unlimited
observer_h
0.0          7.74      11.06
1.7          8.89      12.24
10.0        13.75      17.20
30.0        16.25      19.93

Under eight to nearly twenty square kilometres — a factor of 2.6 — from two parameters. Reporting "the viewshed is 12 km²" without them is reporting a choice, not a measurement.

Note that the range limit costs proportionally less at greater observer heights (30% of the area at ground level, 18% from a mast). A high observer sees more of the near ground too, so a 5 km cap removes a smaller share of what it could see.

Example 3 — a cumulative viewshed from several observers 

def cumulative_viewshed(path, observers, **kwargs):
    """How many observers can see each cell — the standard multi-turbine product."""
    total = None
    for i, xy in enumerate(observers, 1):
        vis, meta, profile = viewshed_from_dem(path, xy, **kwargs)
        total = vis.astype("int16") if total is None else total + vis
        print(f"    observer {i}: {meta['visible_km2']} km²")

    seen_by_any = (total > 0)
    print(f"  seen by at least one: {seen_by_any.mean():.2%}")
    print(f"  seen by all {len(observers)}: {(total == len(observers)).mean():.2%}")
    print(f"  mean visible count where seen: {total[seen_by_any].mean():.2f}")
    return total, profile


turbines = [(-4.0763, 53.0685), (-4.0600, 53.0550), (-4.0900, 53.0800)]
count, profile = cumulative_viewshed("snowdonia_glo30.tif", turbines,
                                     observer_h=100.0, target_h=1.7)
    observer 1: 30.18 km²
    observer 2: 24.61 km²
    observer 3: 26.03 km²
  seen by at least one: 71.44%
  seen by all 3: 18.09%
  mean visible count where seen: 1.62

Note the parameters: observer_h=100 puts the eye at turbine hub height and target_h=1.7 puts the target at a person's eye. That is the reciprocal of the usual formulation, and it is equivalent — visibility is symmetric — but it means one ray cast per turbine rather than one per viewer, which is enormously cheaper.

The individual areas sum to 80.8 km² while the union is 36.9 km² — the three viewsheds overlap heavily. Reporting the sum would overstate the impact by a factor of two.

Explanation

Why the running maximum is the whole algorithm

Along one ray, a cell is hidden if anything closer to the observer subtends a greater vertical angle. So the state you need is a single number — the largest angle seen so far — and it can only increase.

That makes the test angle >= running_max, and np.maximum.accumulate computes the running maximum for a whole ray in one vectorised call. The entire viewshed is then a loop over boundary cells, with everything inside the loop done in NumPy.

For this 252 × 239 DEM that is about 1,000 rays and 0.06 seconds. The naive per-cell approach — checking every intervening cell for every target cell — is O(n²·√n) and takes minutes.

Why this is the R2 algorithm, and what it approximates

Casting rays to boundary cells only is the R2 algorithm. It is fast and it has a known bias: cells near the observer are crossed by many rays and cells far away by few, so distant visibility is slightly under-sampled.

R3, the exact version, casts a ray to every cell. It is exact and far slower — roughly n times more rays.

For most purposes R2 is fine, and its error is a small number of isolated cells near the edge that should have been visible. If your use is legally consequential, use a proper implementation — GDAL's gdal_viewshed, GRASS r.viewshed, or whitebox — which handle this and several other subtleties.

Rays cast from an observer to boundary cells, dense near the observer and sparse at the edges.
R2 casts one ray per boundary cell, so distant cells are sampled less densely than near ones.

Why curvature matters only at range

The earth's surface drops below a tangent plane by approximately d² / (2R):

Distance Drop With refraction
5 km 2.0 m 1.7 m
10 km 7.8 m 6.8 m
30 km 70.6 m 61.4 m
60 km 282 m 246 m

At 5 km the correction is smaller than the DEM's own vertical accuracy, which is why the measured difference over this 8 km extent was 0.05 percentage points. At 60 km it is 246 m — larger than most hills.

Atmospheric refraction bends light downward around the earth, effectively increasing the radius. The standard coefficient of 0.13 recovers about 13% of the drop; it varies with temperature gradient and is much larger over cold water, which is why mirages happen at sea.

Why a DSM and a DTM give different answers

A DTM viewshed says what you could see if the world were bare ground. A DSM viewshed includes trees and buildings as obstructions.

Neither is right for every question. For a turbine assessment, screening by woodland is real and a DSM is closer to the truth — but forestry changes, so a DSM assessment has a shelf life. Planning guidance often specifies a bare-earth viewshed precisely to avoid depending on vegetation that may be felled.

What is never right is being unaware of which one you used. See DEM, DSM and DTM.

Edge cases or notes

  • NoData blocks rays. Fill voids first, or a void casts a shadow that is an artefact.
  • observer_h=0 is almost never intended. Use 1.7 m for a person, or the real structure height.
  • target_h answers a different question from observer_h — what you are looking at, not from.
  • Visibility is symmetric, so a wind-farm viewshed is cheaper computed from each turbine outward than from every viewpoint inward.
  • Curvature is negligible under 10 km and dominant beyond 30 km. The refraction coefficient of 0.13 is a standard-atmosphere value.
  • The DEM must extend beyond the range of interest, or the viewshed is clipped by the grid rather than by terrain.
  • Cumulative viewsheds do not sum. Overlapping areas would be double-counted; take the union.
  • gdal_viewshed, GRASS r.viewshed and whitebox are the production tools. This implementation is for understanding what they do.

FAQ

What observer height should I use?

1.7 m for a standing person, or the actual height of the structure. Zero is rarely what anyone means and understates visibility by about 11% here.

What is the difference between observer height and target height?

Observer height raises the eye; target height raises what you are looking at. For a turbine assessment, target_h is the turbine height and it changes the answer completely.

Does earth curvature matter?

Not under about 10 km — the correction was 0.05 percentage points over an 8 km extent. Beyond 30 km it dominates: the drop is over 70 m.

What is the refraction coefficient?

Atmospheric refraction bends light downward, recovering about 13% of the curvature drop. 0.13 is the standard-atmosphere value; it is much larger over cold water.

Should I use a DSM or a DTM?

A DSM includes trees and buildings as obstructions, which is more realistic but changes as vegetation does. Planning guidance often specifies bare earth for that reason. Record which you used.

How do I combine viewsheds from several points?

Sum the boolean rasters to get a count of how many observers see each cell, then take > 0 for the union. Never add the individual areas — they overlap.

Is this fast enough for real work?

For a few hundred observers on a moderate DEM, yes — about 0.06 s each here. For legally consequential assessments use gdal_viewshed, GRASS r.viewshed or whitebox, which implement the exact algorithm.