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 large margin. From Snowdon's summit, over the Copernicus DEM (step 2 has the code):

observer at the summit, 1,072.8 m
  eye height 0 m                visible 30.24%
  eye height 1.7 m              visible 31.18%
  30 m mast                     visible 38.99%

A 30 m mast sees 25% 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
        ground = np.divide(z - drop - eye, d, out=np.full_like(d, -np.inf), where=d > 0)
        target = np.divide(z + target_h - drop - eye, d, out=np.full_like(d, -np.inf), where=d > 0)
        if max_range:
            target = np.where(d <= max_range, target, -np.inf)

        # the horizon is made by the ground; only the cell being tested is raised by target_h
        horizon = np.maximum.accumulate(np.concatenate([[-np.inf], ground[:-1]]))
        seen = target >= 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 31.18% (18,855 cells)

The angle test, observer and target heights, curvature

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 30.24% (18,290 cells)
  a person        1.7 m: visible 31.18% (18,855 cells)
  a house          10 m: visible 34.74% (21,010 cells)
  a mast           30 m: visible 38.99% (23,582 cells)

The first 1.7 m adds 565 cells — 3% more visible area — because it lifts the eye above the micro-relief of the summit itself; from a summit most of the view is already open, and it takes a mast to add a lot more. From a valley floor the same 1.7 m matters far more.

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

ground_visibility = viewshed(dem, row, col, 27.9, 30.7, observer_h=1.7, target_h=0.0)
turbine_visibility = viewshed(dem, row, col, 27.9, 30.7,
                              observer_h=1.7, target_h=100.0)
print(f"ground visible over {ground_visibility.mean():.2%}; "
      f"the tip of a 100 m turbine visible over {turbine_visibility.mean():.2%}")

ground visible over 31.18%; the tip of a 100 m turbine visible over 58.94%

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 31.18%
  curvature False: visible 31.22%

A difference of 0.04 percentage points — negligible, because this DEM spans about 8 km. The earth drops roughly $d^2/(2R)$ below a tangent plane: 5 m at 8 km, but 80 m at 32 km and 321 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 31.18%
      20 km: visible 31.18%
      10 km: visible 31.18%
       5 km: visible 24.34%

The 5 km limit removes more than a fifth 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.

A recorded viewshed from the summit, and a cumulative one

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, 0.1° × 0.07° around the summit, read from its public bucket; the observer stands on the summit cell. The window crosses the 4°W tile boundary, so two tiles are merged.

import numpy as np
import rasterio
from rasterio.merge import merge
from rasterio.transform import rowcol

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("snowdonia_glo30.tif", "w", **profile) as dst:
    dst.write(mosaic)
dem = mosaic[0].astype("float64")
row, col = rowcol(transform, -4.0763, 53.0685)              # Snowdon summit
print(f"DEM {dem.shape[1]} × {dem.shape[0]}; summit cell ({row}, {col}) at {dem[row, col]:.1f} m")

DEM 240 × 252; summit cell (113, 56) at 1072.8 m

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
        ground = np.divide(z - drop - eye, d, out=np.full_like(d, -np.inf), where=d > 0)
        target = np.divide(z + target_h - drop - eye, d, out=np.full_like(d, -np.inf), where=d > 0)
        if max_range:
            target = np.where(d <= max_range, target, -np.inf)

        # the horizon is made by the ground; only the cell being tested is raised by target_h
        horizon = np.maximum.accumulate(np.concatenate([[-np.inf], ground[:-1]]))
        seen = target >= 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 1072.8 m + 1.7 m, target +0.0 m
  visible 31.18% = 16.14 km² (18,855 cells)

The DEM covers about 51.8 km², so just under a third 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)

  observer 1072.8 m + 0 m, target +0.0 m
  visible 30.24% = 15.66 km² (18,290 cells)
  observer 1072.8 m + 0 m, target +0.0 m
  visible 23.49% = 12.16 km² (14,204 cells)
  observer 1072.8 m + 1.7 m, target +0.0 m
  visible 31.18% = 16.14 km² (18,855 cells)
  observer 1072.8 m + 1.7 m, target +0.0 m
  visible 24.36% = 12.61 km² (14,733 cells)
  observer 1072.8 m + 10 m, target +0.0 m
  visible 34.74% = 17.99 km² (21,010 cells)
  observer 1072.8 m + 10 m, target +0.0 m
  visible 27.73% = 14.36 km² (16,773 cells)
  observer 1072.8 m + 30 m, target +0.0 m
  visible 38.99% = 20.19 km² (23,582 cells)
  observer 1072.8 m + 30 m, target +0.0 m
  visible 31.55% = 16.34 km² (19,082 cells)
range        5 km  unlimited
observer_h                  
0.0         12.16      15.66
1.7         12.61      16.14
10.0        14.36      17.99
30.0        16.34      20.19

From 12.2 to 20.2 square kilometres — a factor of 1.7 — from two parameters, on one summit. Reporting "the viewshed is 16 km²" without them is reporting a choice, not a measurement.

Note that the range limit costs proportionally less at greater observer heights (22% of the area at ground level, 19% 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)]    # hypothetical positions
count, profile = cumulative_viewshed("snowdonia_glo30.tif", turbines,
                                     observer_h=100.0, target_h=1.7)

  observer 1072.8 m + 100.0 m, target +1.7 m
  visible 49.12% = 25.43 km² (29,709 cells)
    observer 1: 25.43 km²
  observer 579.5 m + 100.0 m, target +1.7 m
  visible 21.12% = 10.94 km² (12,776 cells)
    observer 2: 10.94 km²
  observer 590.0 m + 100.0 m, target +1.7 m
  visible 3.57% = 1.85 km² (2,159 cells)
    observer 3: 1.85 km²
  seen by at least one: 54.62%
  seen by all 3: 0.00%
  mean visible count where seen: 1.35

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 nearly equivalent — a sightline works both ways, though the rays are sampled from the other end — but it means one ray cast per turbine rather than one per viewer, which is enormously cheaper.

The individual areas sum to 38.2 km² while the union is about 28.3 km² — the viewsheds overlap. Reporting the sum would overstate the impact by more than a third. The positions are hypothetical, and the summit one dominates: a 100 m structure on Snowdon would be seen from most of this window.

The running maximum along each ray

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.

Why curvature matters only at range

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

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

NoData shadows, symmetry and DEM extent

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

Internal links

• Digital elevation models explained — DSM versus DTM for visibility
• How to extract an elevation profile along a line in Python — the single-ray case
• How to calculate slope and aspect from a DEM in Python — the cell-size handling reused here
• How to make a hillshade from a DEM in Python — for presenting the viewshed over terrain
• Vertical datums explained — why a constant offset does not affect visibility
• The raster data model explained — writing the boolean output
• How to calculate zonal statistics in Python — how much of each parish can see the site
• How to clip a raster to a polygon in Python — restricting the viewshed to a study area
• How to test whether a window can see a landmark — line of sight through buildings

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; from Snowdon's summit it understates visibility by 3%, and from lower ground, where the eye is closer to the obstructions, by more.

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