Source code: Viewshed analysis in Python: line of sight and a wind farm ZTV with GDAL and numpy

"""Viewshed analysis in Python: a zone of theoretical visibility (ZTV) for Scout Moor wind farm,
Greater Manchester / Lancashire, with GDAL's ViewshedGenerate and a from-scratch numpy line-of-sight
viewshed. Covers earth curvature and refraction, a cumulative viewshed from 26 turbines, line-of-sight
profiles, and how much a LiDAR DSM (trees and buildings) shrinks a bare-earth DTM viewshed.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/viewshed-analysis-python/
Data: Copernicus DEM GLO-30 (ESA / Copernicus, free licence, via the Microsoft Planetary Computer);
Environment Agency LIDAR Composite DTM and first-return DSM 1 m (Open Government Licence v3.0);
turbine positions (c) OpenStreetMap contributors (ODbL).
"""
import json
import sys
import time
import warnings
from importlib.metadata import version
from pathlib import Path

import numpy as np
import requests
import pystac_client
import planetary_computer
from osgeo import gdal, osr
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm, ListedColormap

gdal.UseExceptions()
warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)

TIP_H, HUB_H = 100.0, 60.0       # Nordex N80 at Scout Moor: 60 m hub + 40 m blade = 100 m to blade tip
EYE_H = 2.0                      # viewer eye height used in UK landscape guidance (about 1.5-2 m)
RADIUS = 30_000                  # ZTV radius for turbines of 86-100 m to blade tip (NatureScot guidance)
CELL = 30.0                      # GLO-30 is 1 arc-second (about 30 m); we warp it to a 30 m BNG grid
CURV = 0.85714                   # 1 - k with refraction coefficient k = 1/7: the GDAL default for visible light
R_EARTH = 6_371_000.0
SPECTRAL = plt.get_cmap("Spectral_r")
PLACES = {                       # name: (lat, lon) of town centres and viewpoints
    "Manchester Town Hall": (53.4794, -2.2446), "Bury": (53.5933, -2.2966), "Rochdale": (53.6163, -2.1561),
    "Bolton": (53.5781, -2.4290), "Oldham": (53.5409, -2.1114), "Ramsbottom": (53.6476, -2.3166),
    "Edenfield": (53.6705, -2.3062), "Rawtenstall": (53.7000, -2.2890), "Burnley": (53.7893, -2.2405),
    "Blackburn": (53.7486, -2.4823), "Peel Monument": (53.6434, -2.3380), "Winter Hill": (53.6298, -2.5159),
}


def section(title):
    print(f"\n=== {title} ===")


section("versions")
print(f"python {sys.version.split()[0]} | GDAL {gdal.__version__} | numpy {np.__version__} | "
      f"matplotlib {matplotlib.__version__} | pystac-client {version('pystac-client')}")

BNG = osr.SpatialReference(); BNG.ImportFromEPSG(27700)
WGS84 = osr.SpatialReference(); WGS84.ImportFromEPSG(4326)
WGS84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)      # lon, lat order
TO_BNG = osr.CoordinateTransformation(WGS84, BNG)


def to_bng(lat, lon):
    x, y, _ = TO_BNG.TransformPoint(lon, lat)
    return x, y


# ---------------------------------------------------------------- Step 1: turbine positions from OpenStreetMap
section("Step 1: Scout Moor turbines from OpenStreetMap")
OVERPASS = ["https://overpass-api.de/api/interpreter", "https://overpass.private.coffee/api/interpreter",
            "https://maps.mail.ru/osm/tools/overpass/api/interpreter"]
QUERY = """[out:json][timeout:60];
node["power"="generator"]["generator:source"="wind"]["generator:output:electricity"="2.5 MW"]
  (53.62,-2.32,53.68,-2.18);
out;"""


def overpass(query, cache):
    if cache.exists():
        return json.loads(cache.read_text())
    for attempt in range(4):
        for url in OVERPASS:
            try:
                r = requests.post(url, data={"data": query}, timeout=90,
                                  headers={"User-Agent": "technicalgis-tutorial"})
                r.raise_for_status()
                els = r.json()["elements"]
                cache.write_text(json.dumps(els))
                return els
            except (requests.RequestException, ValueError):
                continue
        time.sleep(15 * (attempt + 1))
    raise RuntimeError("all Overpass mirrors failed")


turbines = overpass(QUERY, DATA / "scout_moor_turbines.json")
tx, ty = np.array([to_bng(e["lat"], e["lon"]) for e in turbines]).T
print(f"turbines: {len(turbines)} | E {tx.min():.0f}..{tx.max():.0f} | N {ty.min():.0f}..{ty.max():.0f}")
cx, cy = round(tx.mean(), -1), round(ty.mean(), -1)
spacing = np.sort(np.hypot(tx[:, None] - tx, ty[:, None] - ty), axis=1)[:, 1]
print(f"wind farm centre E {cx:.0f} N {cy:.0f} | nearest-neighbour spacing {spacing.min():.0f}..{spacing.max():.0f} m")

# ---------------------------------------------------------------- Step 2: Copernicus GLO-30 on a 30 m BNG grid
section("Step 2: Copernicus DEM GLO-30 from the Planetary Computer")
HALF = RADIUS + 1020                       # a little margin around the 30 km circle
BOUNDS = (cx - HALF, cy - HALF, cx + HALF, cy + HALF)
dem_path = DATA / "glo30_bng_30m.tif"
if not dem_path.exists():
    catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                                        modifier=planetary_computer.sign_inplace)
    items = list(catalog.search(collections=["cop-dem-glo-30"], bbox=[-2.79, 53.33, -1.69, 53.98]).items())
    print("tiles:", ", ".join(sorted(it.id for it in items)))
    gdal.Warp(str(dem_path), ["/vsicurl/" + it.assets["data"].href for it in items], dstSRS="EPSG:27700",
              outputBounds=BOUNDS, xRes=CELL, yRes=CELL, resampleAlg="bilinear", outputType=gdal.GDT_Float32,
              creationOptions=["COMPRESS=DEFLATE", "TILED=YES"])
ds = gdal.Open(str(dem_path))
dem = ds.GetRasterBand(1).ReadAsArray()
GT = ds.GetGeoTransform()
print(f"grid {ds.RasterXSize} x {ds.RasterYSize} cells of {GT[1]:.0f} m | elevation {dem.min():.1f}..{dem.max():.1f} m")


def rowcol(x, y, gt=GT):
    return int((gt[3] - y) // -gt[5]), int((x - gt[0]) // gt[1])


tr, tc = np.array([rowcol(x, y) for x, y in zip(tx, ty)]).T
print(f"turbine base elevation (GLO-30): {dem[tr, tc].min():.0f}..{dem[tr, tc].max():.0f} m, "
      f"tips at {dem[tr, tc].min() + TIP_H:.0f}..{dem[tr, tc].max() + TIP_H:.0f} m")
yy, xx = np.mgrid[0:dem.shape[0], 0:dem.shape[1]]
dist = np.hypot(GT[0] + (xx + 0.5) * GT[1] - cx, GT[3] + (yy + 0.5) * GT[5] - cy)
in_circle = dist <= RADIUS
CELL_KM2 = CELL * CELL / 1e6
print(f"30 km study circle: {in_circle.sum() * CELL_KM2:,.0f} km2")


def hillshade(z, cell, az=315, alt=45):
    gy, gx = np.gradient(z, cell)
    slope = np.arctan(np.hypot(gx, gy))
    aspect = np.arctan2(-gx, gy)
    a, b = np.radians(360 - az + 90), np.radians(alt)
    return np.clip(np.sin(b) * np.cos(slope) + np.cos(b) * np.sin(slope) * np.cos(a - aspect), 0, 1)


hs = hillshade(dem, CELL)
EXTENT = [(GT[0] - cx) / 1000, (GT[0] + dem.shape[1] * GT[1] - cx) / 1000,
          (GT[3] + dem.shape[0] * GT[5] - cy) / 1000, (GT[3] - cy) / 1000]
place_xy = {k: to_bng(*v) for k, v in PLACES.items()}

fig, ax = plt.subplots(figsize=(7.2, 6.4))
im = ax.imshow(dem, cmap="Spectral_r", extent=EXTENT, vmin=0, vmax=560)
ax.imshow(hs, cmap="gray", alpha=0.35, extent=EXTENT)
ax.add_patch(plt.Circle((0, 0), RADIUS / 1000, fill=False, ls="--", lw=1, color="k"))
ax.plot((tx - cx) / 1000, (ty - cy) / 1000, "^", ms=4, mfc="#d53e4f", mec="k", mew=0.4, label="Scout Moor turbines (26)")
for name, (x, y) in place_xy.items():
    left = name in ("Ramsbottom", "Peel Monument")    # nudge labels that would sit on top of the turbines
    ax.plot((x - cx) / 1000, (y - cy) / 1000, "o", ms=3, color="k")
    ax.annotate(name, ((x - cx) / 1000, (y - cy) / 1000), xytext=(3, -9) if name == "Ramsbottom" else
                ((-3, 3) if left else (3, 3)), textcoords="offset points", fontsize=7,
                ha="right" if name == "Peel Monument" else "left")
ax.set_xlabel("km east of the wind farm centre"); ax.set_ylabel("km north")
ax.set_title("Copernicus GLO-30 around Scout Moor wind farm, with the 30 km ZTV radius", loc="left")
ax.legend(loc="lower left", fontsize=7)
fig.colorbar(im, ax=ax, shrink=0.8, label="elevation (m, EGM2008)")
fig.savefig(OUT / "01-study-area-dem.png"); plt.close(fig)

# ---------------------------------------------------------------- Step 3: one viewshed with GDAL
section("Step 3: single-turbine viewshed with gdal.ViewshedGenerate")
band = ds.GetRasterBand(1)


def gdal_viewshed(band, x, y, obs_h, tgt_h, curv=CURV, max_dist=RADIUS, height_mode=gdal.GVOT_NORMAL,
                  shape=dem.shape, gt=GT):
    """Run GDAL's viewshed in memory and paste the result (cropped to max_dist) onto the full grid."""
    oor = 0 if height_mode == gdal.GVOT_NORMAL else 1e9     # out of range: never visible
    out = gdal.ViewshedGenerate(band, "MEM", "", [], x, y, obs_h, tgt_h, 1, 0, oor, -1, curv,
                                gdal.GVM_Edge, max_dist, heightMode=height_mode)
    a = out.GetRasterBand(1).ReadAsArray()
    ogt = out.GetGeoTransform()
    r0, c0 = int(round((gt[3] - ogt[3]) / -gt[5])), int(round((ogt[0] - gt[0]) / gt[1]))
    fill = 1e9 if height_mode != gdal.GVOT_NORMAL else 0
    full = np.full(shape, fill, dtype="float32")
    full[r0:r0 + a.shape[0], c0:c0 + a.shape[1]] = a
    return full


T = int(np.argmin(np.hypot(tx - cx, ty - cy)))          # the turbine closest to the centre of the farm
t0 = time.perf_counter()
vs_gdal = gdal_viewshed(band, tx[T], ty[T], TIP_H, EYE_H) == 1
t_gdal = time.perf_counter() - t0
vis_km2 = (vs_gdal & in_circle).sum() * CELL_KM2
print(f"turbine T{T + 1} at E {tx[T]:.0f} N {ty[T]:.0f}, base {dem[tr[T], tc[T]]:.0f} m, observer = blade tip {TIP_H:.0f} m")
print(f"GDAL viewshed: {t_gdal:.2f} s | tip visible from {vis_km2:,.0f} km2 = "
      f"{100 * vis_km2 / (in_circle.sum() * CELL_KM2):.1f}% of the 30 km circle")

# ---------------------------------------------------------------- Step 4: the same viewshed in numpy
section("Step 4: a numpy line-of-sight viewshed (ray casting)")


def numpy_viewshed(z, r0, c0, obs_h, tgt_h, cell, curv=CURV, max_dist=RADIUS):
    """Cast one ray to every edge cell of a square around the observer. Along each ray a cell is visible
    if the slope from the eye to the target above it is at least the steepest terrain slope seen so far."""
    n = int(max_dist / cell) + 1
    edge = np.arange(-n, n + 1)
    ends = np.concatenate([np.c_[np.full_like(edge, -n), edge], np.c_[np.full_like(edge, n), edge],
                           np.c_[edge, np.full_like(edge, -n)], np.c_[edge, np.full_like(edge, n)]])
    k = np.arange(1, n + 1) / n                              # n steps: at most one cell per step in each axis
    rr = np.rint(r0 + ends[:, :1] * k).astype(int)           # rays x steps
    cc = np.rint(c0 + ends[:, 1:] * k).astype(int)
    ok = (rr >= 0) & (rr < z.shape[0]) & (cc >= 0) & (cc < z.shape[1])
    rr, cc = np.clip(rr, 0, z.shape[0] - 1), np.clip(cc, 0, z.shape[1] - 1)
    d = np.hypot(rr - r0, cc - c0) * cell
    drop = curv * d ** 2 / (2 * R_EARTH)                     # earth curvature, reduced by refraction
    z_eye = z[r0, c0] + obs_h
    with np.errstate(divide="ignore", invalid="ignore"):
        s_ground = (z[rr, cc] - drop - z_eye) / d
        s_target = (z[rr, cc] - drop + tgt_h - z_eye) / d
    s_ground[d == 0] = -np.inf
    horizon = np.maximum.accumulate(s_ground, axis=1)
    horizon = np.concatenate([np.full((len(ends), 1), -np.inf), horizon[:, :-1]], axis=1)
    visible = (s_target >= horizon) & ok & (d <= max_dist)
    flat = (rr * z.shape[1] + cc)[ok]
    hits = np.bincount(flat, weights=visible[ok], minlength=z.size)   # any ray that sees the cell wins
    out = (hits > 0).reshape(z.shape)
    out[r0, c0] = True
    return out


t0 = time.perf_counter()
vs_np = numpy_viewshed(dem, tr[T], tc[T], TIP_H, EYE_H, CELL)
t_np = time.perf_counter() - t0
both = in_circle
agree = (vs_np == vs_gdal)[both].mean()
pe = vs_np[both].mean() * vs_gdal[both].mean() + (~vs_np[both]).mean() * (~vs_gdal[both]).mean()
kappa = (agree - pe) / (1 - pe)
print(f"numpy viewshed: {t_np:.2f} s | tip visible from {(vs_np & both).sum() * CELL_KM2:,.0f} km2")
print(f"agreement with GDAL: {100 * agree:.2f}% of cells, Cohen's kappa {kappa:.3f} | "
      f"numpy-only {((vs_np & ~vs_gdal) & both).sum() * CELL_KM2:.1f} km2, "
      f"GDAL-only {((vs_gdal & ~vs_np) & both).sum() * CELL_KM2:.1f} km2")

zoom = 12  # km half-width for the comparison maps
fig, axs = plt.subplots(1, 3, figsize=(12, 4.5), sharex=True, sharey=True, layout="constrained")
diff = np.zeros(dem.shape, "int8"); diff[vs_gdal & vs_np] = 1; diff[vs_gdal & ~vs_np] = 2; diff[~vs_gdal & vs_np] = 3
panels = [(vs_gdal, f"gdal.ViewshedGenerate ({t_gdal:.2f} s)"), (vs_np, f"numpy ray casting ({t_np:.2f} s)")]
for ax, (v, title) in zip(axs[:2], panels):
    ax.imshow(hs, cmap="gray", extent=EXTENT, vmin=0, vmax=1)
    ax.imshow(np.ma.masked_where(~v, v), cmap=ListedColormap(["#f46d43"]), alpha=0.7, extent=EXTENT)
    ax.set_title(title, loc="left")
cm = ListedColormap(["#3288bd", "#d53e4f", "#fdae61"])
axs[2].imshow(hs, cmap="gray", extent=EXTENT, vmin=0, vmax=1)
axs[2].imshow(np.ma.masked_where(diff == 0, diff), cmap=cm, norm=BoundaryNorm([0.5, 1.5, 2.5, 3.5], 3),
              alpha=0.75, extent=EXTENT, interpolation="nearest")
axs[2].set_title(f"agreement {100 * agree:.1f}% (kappa {kappa:.2f})", loc="left")
for lab, col in [("both visible", "#3288bd"), ("GDAL only", "#d53e4f"), ("numpy only", "#fdae61")]:
    axs[2].plot([], [], "s", color=col, label=lab)
axs[2].legend(loc="lower left", fontsize=7)
for ax in axs:
    ax.plot((tx[T] - cx) / 1000, (ty[T] - cy) / 1000, "^", ms=6, mfc="k", mec="w")
    ax.set_xlim(-zoom, zoom); ax.set_ylim(-zoom, zoom); ax.set_xlabel("km east")
axs[0].set_ylabel("km north")
fig.suptitle(f"Where the blade tip of one turbine (T{T + 1}, {TIP_H:.0f} m) can be seen from 2 m eye height",
             x=0.01, ha="left")
fig.savefig(OUT / "02-viewshed-gdal-vs-numpy.png"); plt.close(fig)

# ---------------------------------------------------------------- Step 5: earth curvature and refraction
section("Step 5: earth curvature and refraction")
for d_km in (5, 10, 20, 30):
    d = d_km * 1000
    print(f"{d_km:>2} km: curvature drop {d ** 2 / (2 * R_EARTH):6.1f} m | with refraction (k=1/7) "
          f"{CURV * d ** 2 / (2 * R_EARTH):6.1f} m")
coeffs = {"flat earth (0)": 0.0, "curvature only (1)": 1.0, "curvature + refraction (0.857)": CURV}
bands = np.arange(0, RADIUS + 1, 2500)
band_share = {}
for label, c in coeffs.items():
    v = gdal_viewshed(band, tx[T], ty[T], TIP_H, EYE_H, curv=c) == 1
    km2 = (v & in_circle).sum() * CELL_KM2
    band_share[label] = [v[(dist >= a) & (dist < b)].mean() * 100 for a, b in zip(bands[:-1], bands[1:])]
    far = v[(dist >= 20_000) & (dist < RADIUS)].mean() * 100
    print(f"{label:<32} visible {km2:7.1f} km2 | {far:4.1f}% of the 20-30 km ring")

fig, axs = plt.subplots(1, 2, figsize=(11, 3.8), gridspec_kw={"width_ratios": [1.5, 1]})
cols = ["#3288bd", "#d53e4f", "#66c2a5"]
mid = (bands[:-1] + bands[1:]) / 2000
for (label, share), col in zip(band_share.items(), cols):
    axs[0].plot(mid, share, "o-", color=col, label=label, ms=3)
axs[0].set_xlabel("distance from turbine (km)"); axs[0].set_ylabel("share of ring where the tip is visible (%)")
axs[0].set_title("Visibility by distance band for three curvature settings", loc="left")
axs[0].legend(fontsize=7); axs[0].grid(alpha=0.3)
dd = np.linspace(0, 40, 200)
axs[1].plot(dd, (dd * 1000) ** 2 / (2 * R_EARTH), color="#d53e4f", label="geometric curvature")
axs[1].plot(dd, CURV * (dd * 1000) ** 2 / (2 * R_EARTH), color="#66c2a5", label="with refraction, k = 1/7")
axs[1].set_xlabel("distance (km)"); axs[1].set_ylabel("apparent drop of the ground (m)")
axs[1].set_title("How far the earth curves away", loc="left"); axs[1].legend(fontsize=7); axs[1].grid(alpha=0.3)
fig.savefig(OUT / "03-curvature-refraction.png"); plt.close(fig)

# ---------------------------------------------------------------- Step 6: cumulative ZTV from 26 turbines
section("Step 6: cumulative viewshed (ZTV) from all 26 turbines")
t0 = time.perf_counter()
n_tip = np.zeros(dem.shape, "int16"); n_hub = np.zeros(dem.shape, "int16")
for x, y in zip(tx, ty):
    n_tip += gdal_viewshed(band, x, y, TIP_H, EYE_H) == 1
    n_hub += gdal_viewshed(band, x, y, HUB_H, EYE_H) == 1
print(f"{2 * len(tx)} viewsheds in {time.perf_counter() - t0:.1f} s")
area = in_circle.sum() * CELL_KM2
for label, cond in [("at least one blade tip", n_tip >= 1), ("at least one hub", n_hub >= 1),
                    ("blade tips only, no hub", (n_tip >= 1) & (n_hub == 0)), ("all 26 tips", n_tip == 26),
                    ("no part of any turbine", n_tip == 0)]:
    km2 = (cond & in_circle).sum() * CELL_KM2
    print(f"{label:<26} {km2:7.0f} km2  {100 * km2 / area:5.1f}%")
print("\nturbines visible at named places (2 m eye height, GLO-30 surface):")
print(f"{'place':<22}{'km':>6}{'tips':>6}{'hubs':>6}")
for name, (x, y) in place_xy.items():
    r, c = rowcol(x, y)
    print(f"{name:<22}{np.hypot(x - cx, y - cy) / 1000:6.1f}{n_tip[r, c]:6d}{n_hub[r, c]:6d}")

fig, axs = plt.subplots(1, 2, figsize=(12, 5.4), sharex=True, sharey=True, layout="constrained")
cmap = ListedColormap(SPECTRAL(np.linspace(0.1, 1, 6)))
norm = BoundaryNorm([0.5, 1.5, 5.5, 10.5, 15.5, 20.5, 26.5], cmap.N)
for ax in axs:
    ax.imshow(hs, cmap="gray", extent=EXTENT, vmin=0, vmax=1)
    ax.add_patch(plt.Circle((0, 0), RADIUS / 1000, fill=False, ls="--", lw=0.8, color="k"))
im = axs[0].imshow(np.ma.masked_where((n_tip == 0) | ~in_circle, n_tip), cmap=cmap, norm=norm, extent=EXTENT,
                   alpha=0.85, interpolation="nearest")
cb = fig.colorbar(im, ax=axs[0], shrink=0.75, ticks=[1, 3.5, 8, 13, 18, 23.5])
cb.ax.set_yticklabels(["1", "2-5", "6-10", "11-15", "16-20", "21-26"]); cb.set_label("number of blade tips visible")
axs[0].set_title("Cumulative ZTV: how many of the 26 turbines you can see", loc="left")
cls = np.zeros(dem.shape, "int8"); cls[n_tip >= 1] = 1; cls[n_hub >= 1] = 2
cls = np.ma.masked_where((cls == 0) | ~in_circle, cls)
axs[1].imshow(cls, cmap=ListedColormap(["#fee08b", "#d53e4f"]), norm=BoundaryNorm([0.5, 1.5, 2.5], 2),
              extent=EXTENT, alpha=0.85, interpolation="nearest")
axs[1].plot([], [], "s", color="#fee08b", label="blade tips only")
axs[1].plot([], [], "s", color="#d53e4f", label="hub(s) visible")
axs[1].legend(loc="lower left", fontsize=7)
axs[1].set_title("Blade tips only versus hubs visible", loc="left")
for ax in axs:
    ax.plot((tx - cx) / 1000, (ty - cy) / 1000, "^", ms=2.5, mfc="k", mec="k")
    for name, (x, y) in place_xy.items():
        if np.hypot(x - cx, y - cy) < 6000:      # too close to the farm to label legibly at this scale
            continue
        ax.plot((x - cx) / 1000, (y - cy) / 1000, "o", ms=2.5, color="k")
        ax.annotate(name, ((x - cx) / 1000, (y - cy) / 1000), xytext=(3, 2), textcoords="offset points", fontsize=6)
    ax.set_xlabel("km east")
axs[0].set_ylabel("km north")
fig.savefig(OUT / "04-cumulative-ztv.png"); plt.close(fig)
for name, arr in [("ztv_tips_count", n_tip), ("ztv_hubs_count", n_hub)]:
    drv = gdal.GetDriverByName("GTiff").Create(str(OUT / f"{name}.tif"), dem.shape[1], dem.shape[0], 1,
                                                 gdal.GDT_Int16, ["COMPRESS=DEFLATE"])
    drv.SetGeoTransform(GT); drv.SetProjection(BNG.ExportToWkt()); drv.GetRasterBand(1).WriteArray(arr); drv = None

# ---------------------------------------------------------------- Step 7: line-of-sight profiles
section("Step 7: line-of-sight profiles")


def los_profile(x0, y0, x1, y1, obs_h, tgt_h, z=dem, gt=GT, curv=CURV, step=CELL / 2):
    """Sample the DEM along a straight line, lower it for curvature and test the sightline."""
    L = np.hypot(x1 - x0, y1 - y0)
    d = np.arange(0, L + step, step); d[-1] = L
    xs, ys = x0 + (x1 - x0) * d / L, y0 + (y1 - y0) * d / L
    fc = (xs - gt[0]) / gt[1] - 0.5; fr = (ys - gt[3]) / gt[5] - 0.5            # bilinear sampling
    c0, r0 = np.floor(fc).astype(int), np.floor(fr).astype(int); fx, fy = fc - c0, fr - r0
    g = (z[r0, c0] * (1 - fx) * (1 - fy) + z[r0, c0 + 1] * fx * (1 - fy) + z[r0 + 1, c0] * (1 - fx) * fy
         + z[r0 + 1, c0 + 1] * fx * fy)
    g_eff = g - curv * d ** 2 / (2 * R_EARTH)
    z0, z1 = g_eff[0] + obs_h, g_eff[-1] + tgt_h
    sight = z0 + (z1 - z0) * d / L
    clearance = sight[1:-1] - g_eff[1:-1]
    worst = int(np.argmin(clearance)) + 1
    return d, g, g_eff, sight, clearance.min() >= 0, worst


fig, axs = plt.subplots(2, 1, figsize=(11, 6.4))
for ax, name in zip(axs, ["Manchester Town Hall", "Peel Monument"]):
    x1, y1 = place_xy[name]
    T2 = int(np.argmin(np.hypot(tx - x1, ty - y1)))              # nearest turbine to the viewer
    d, g, g_eff, sight, ok, worst = los_profile(tx[T2], ty[T2], x1, y1, TIP_H, EYE_H)
    margin = (sight - g_eff)[worst]
    print(f"T{T2 + 1} tip -> {name}: {d[-1] / 1000:.1f} km | {'VISIBLE' if ok else 'BLOCKED'} | "
          f"tightest point {d[worst] / 1000:.2f} km from the turbine, sightline {margin:+.1f} m above the "
          f"curvature-corrected ground | curvature + refraction at the viewer {CURV * d[-1] ** 2 / (2 * R_EARTH):.1f} m")
    km = d / 1000
    col = "#1a9850" if ok else "#d53e4f"
    ax.fill_between(km, g_eff, g_eff.min() - 20, color="#abdda4", alpha=0.6, label="GLO-30 surface, lowered for curvature")
    ax.plot(km, g, color="#3288bd", lw=0.8, ls="--", label="GLO-30 surface on a flat earth")
    ax.plot(km, sight, color=col, lw=1.4, label="sightline (clear)" if ok else "sightline (blocked)")
    ax.plot([0, 0], [g_eff[0], g_eff[0] + TIP_H], color="k", lw=2); ax.plot(0, g_eff[0] + TIP_H, "^", color="k")
    ax.plot(km[-1], g_eff[-1] + EYE_H, "o", color="k")
    ax.plot(km[worst], g_eff[worst], "x", color=col, ms=8, mew=2)
    ax.set_xlim(-0.2, km[-1] + 0.2); ax.set_ylim(g_eff.min() - 20, max(g.max(), sight.max()) + 40)
    ax.set_title(f"Turbine T{T2 + 1} blade tip to {name} ({km[-1]:.1f} km): "
                 f"{'visible' if ok else 'hidden'}", loc="left")
    ax.set_ylabel("elevation (m)"); ax.grid(alpha=0.3); ax.legend(fontsize=7, loc="upper right")
axs[1].set_xlabel("distance from the turbine (km)")
fig.savefig(OUT / "05-line-of-sight-profiles.png"); plt.close(fig)

# ---------------------------------------------------------------- Step 8: DTM vs DSM with Environment Agency LiDAR
section("Step 8: bare-earth DTM versus surface DSM (Environment Agency LiDAR, 2 m)")
EA_WCS = {
    "dtm": ("lidar-composite-digital-terrain-model-dtm-1m",
            "13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m"),
    "dsm": ("lidar-composite-digital-surface-model-first-return-dsm-1m",
            "df4e3ec3-315e-48aa-aaaf-b5ae74d7b2bb__Lidar_Composite_Elevation_FZ_DSM_1m"),
}
LOCAL = (378000, 414000, 386000, 422000)    # 8 km x 8 km: Ramsbottom, Edenfield and the western half of the farm
RES = 2.0


def ea_tile(product, x0, y0):
    path = DATA / "ea" / f"{product}_{x0}_{y0}.tif"
    if not path.exists():
        path.parent.mkdir(exist_ok=True)
        service, cov = EA_WCS[product]
        url = (f"https://environment.data.gov.uk/spatialdata/{service}/wcs?service=WCS&version=2.0.1"
               f"&request=GetCoverage&CoverageId={cov}&format=image/tiff"
               f"&subset=E({x0},{x0 + 1000})&subset=N({y0},{y0 + 1000})")
        for attempt in range(4):
            try:
                r = requests.get(url, timeout=180); r.raise_for_status(); break
            except requests.RequestException:
                if attempt == 3:
                    raise
                time.sleep(5)
        path.write_bytes(r.content)
    return str(path)


lidar = {}
for product in ("dtm", "dsm"):
    out = DATA / f"ea_{product}_2m.tif"
    if not out.exists():
        tiles = [ea_tile(product, x0, y0) for x0 in range(LOCAL[0], LOCAL[2], 1000)
                 for y0 in range(LOCAL[1], LOCAL[3], 1000)]
        gdal.Warp(str(out), tiles, outputBounds=LOCAL, xRes=RES, yRes=RES, resampleAlg="average",
                  srcNodata=-3.4028234663852886e+38, dstNodata=np.nan, outputType=gdal.GDT_Float32)
    src = gdal.Open(str(out))              # keep the dataset alive while its band is read
    lidar[product], LGT = src.GetRasterBand(1).ReadAsArray(), src.GetGeoTransform()
    src = None
dtm, dsm = lidar["dtm"], lidar["dsm"]
print(f"grid {dtm.shape[1]} x {dtm.shape[0]} cells of {RES:.0f} m | no data: DTM {np.isnan(dtm).mean() * 100:.2f}%, "
      f"DSM {np.isnan(dsm).mean() * 100:.2f}%")
dsm = np.where(np.isnan(dsm), dtm, dsm); dtm = np.where(np.isnan(dtm), dsm, dtm)
dsm = np.fmax(dsm, dtm)
ndsm = dsm - dtm
local_turb = [(x, y) for x, y in zip(tx, ty) if LOCAL[0] < x < LOCAL[2] and LOCAL[1] < y < LOCAL[3]]
print(f"turbines inside the LiDAR window: {len(local_turb)} of {len(tx)}")
# The DSM contains the towers themselves: put the observer back on the ground at each turbine
LY, LX = np.mgrid[0:dtm.shape[0], 0:dtm.shape[1]]
tower = []
for x, y in local_turb:
    r, c = int((LGT[3] - y) // RES), int((x - LGT[0]) // RES)
    near = (LY - r) ** 2 + (LX - c) ** 2 <= 6 ** 2
    tower.append(ndsm[near].max())
    dsm[near] = dtm[near]
print(f"DSM - DTM within 12 m of each turbine (the tower in the DSM): median {np.median(tower):.1f} m, "
      f"max {np.max(tower):.1f} m; reset to the DTM")
ndsm = dsm - dtm
open_ground = ndsm < 1.0
print(f"cells more than 2.5 m above ground (trees, buildings): {(ndsm > 2.5).mean() * 100:.1f}%")


def mem_band(a, gt):
    d = gdal.GetDriverByName("MEM").Create("", a.shape[1], a.shape[0], 1, gdal.GDT_Float32)
    d.SetGeoTransform(gt); d.SetProjection(BNG.ExportToWkt()); d.GetRasterBand(1).WriteArray(a)
    return d


dtm_ds, dsm_ds = mem_band(dtm, LGT), mem_band(dsm, LGT)
n_dtm = np.zeros(dtm.shape, "int16"); n_dsm = np.zeros(dtm.shape, "int16")
t0 = time.perf_counter()
for x, y in local_turb:
    n_dtm += gdal_viewshed(dtm_ds.GetRasterBand(1), x, y, TIP_H, EYE_H, max_dist=15_000,
                           shape=dtm.shape, gt=LGT) == 1
    # lowest elevation a target in each cell must reach to be seen over the DSM; a ground-level eye is at DTM + 2 m
    need = gdal_viewshed(dsm_ds.GetRasterBand(1), x, y, TIP_H, 0, max_dist=15_000, shape=dtm.shape, gt=LGT,
                         height_mode=gdal.GVOT_MIN_TARGET_HEIGHT_FROM_DEM)
    n_dsm += np.nan_to_num(need <= dtm + EYE_H, nan=0).astype("int16")
print(f"{2 * len(local_turb)} LiDAR viewsheds in {time.perf_counter() - t0:.1f} s")
for label, m in [("all cells", np.ones(dtm.shape, bool)), ("open ground only", open_ground)]:
    a_dtm, a_dsm = (n_dtm >= 1)[m].mean() * 100, (n_dsm >= 1)[m].mean() * 100
    print(f"{label:<17}: >=1 tip visible on DTM {a_dtm:5.1f}% | on DSM {a_dsm:5.1f}% | "
          f"mean tips visible {n_dtm[m].mean():4.1f} vs {n_dsm[m].mean():4.1f}")
for name in ("Edenfield", "Ramsbottom"):
    x, y = place_xy[name]
    r, c = int((LGT[3] - y) // RES), int((x - LGT[0]) // RES)
    win = (slice(max(r - 250, 0), r + 250), slice(max(c - 250, 0), c + 250))     # 1 km x 1 km around the centre
    print(f"{name:<11} 1 km2 around ({x:.0f}, {y:.0f}): >=1 tip visible on DTM {(n_dtm[win] >= 1).mean() * 100:5.1f}% | "
          f"on DSM (at ground level) {(n_dsm[win] >= 1).mean() * 100:5.1f}%")

LEXT = [(LOCAL[0] - cx) / 1000, (LOCAL[2] - cx) / 1000, (LOCAL[1] - cy) / 1000, (LOCAL[3] - cy) / 1000]
lhs = hillshade(dsm, RES)
fig, axs = plt.subplots(1, 3, figsize=(13, 4.6), sharex=True, sharey=True, layout="constrained")
axs[0].imshow(lhs, cmap="gray", extent=LEXT, vmin=0, vmax=1)
im0 = axs[0].imshow(np.ma.masked_less(ndsm, 2.5), cmap="Spectral_r", vmin=2.5, vmax=25, extent=LEXT, alpha=0.8)
fig.colorbar(im0, ax=axs[0], shrink=0.7, label="DSM - DTM (m)")
axs[0].set_title("LiDAR surface objects over 2.5 m", loc="left")
cmap6 = ListedColormap(SPECTRAL(np.linspace(0.1, 1, 6)))
nrm = BoundaryNorm([0.5, 1.5, 3.5, 6.5, 9.5, 12.5, 30], cmap6.N)
for ax, n, title in [(axs[1], n_dtm, "Bare-earth DTM"), (axs[2], n_dsm, "DSM, viewer at ground level")]:
    ax.imshow(lhs, cmap="gray", extent=LEXT, vmin=0, vmax=1)
    im = ax.imshow(np.ma.masked_equal(n, 0), cmap=cmap6, norm=nrm, extent=LEXT, alpha=0.85, interpolation="nearest")
    ax.set_title(f"{title}: {(n >= 1).mean() * 100:.0f}% see a tip", loc="left")
cb = fig.colorbar(im, ax=axs[2], shrink=0.7, ticks=[1, 2.5, 5, 8, 11, 14])
cb.ax.set_yticklabels(["1", "2-3", "4-6", "7-9", "10-12", "13+"]); cb.set_label("blade tips visible")
for ax in axs:
    ax.plot([(x - cx) / 1000 for x, _ in local_turb], [(y - cy) / 1000 for _, y in local_turb], "^",
            ms=3, mfc="k", mec="w", mew=0.3)
    for name in ("Edenfield", "Ramsbottom", "Peel Monument"):
        x, y = place_xy[name]
        ax.annotate(name, ((x - cx) / 1000, (y - cy) / 1000), fontsize=7, ha="center",
                    bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.75))
    ax.set_xlabel("km east")
axs[0].set_ylabel("km north")
fig.savefig(OUT / "06-dtm-vs-dsm-lidar.png"); plt.close(fig)
print("\nfigures written to", OUT.resolve())