Voronoi and Delaunay: the rainfall over a catchment from twelve gauges

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patheffects import withStroke

plt.rcParams.update({
    "figure.figsize": (7, 4.8), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"
RAIN = LinearSegmentedColormap.from_list("rain", ["#ffffff", SECOND])   # one hue, light to dark

# The catchment, the gauges, and the rain are synthetic, so that the answer can be
# checked at the end. Pretend you have not seen this cell.
catchment = np.array([                                    # outline in local plane coordinates, km
    (0.5, 2.5), (4.5, 0.3), (9.0, 0.8), (13.5, 0.0), (18.0, 1.5), (22.0, 3.0), (24.0, 7.0),
    (23.5, 12.0), (21.0, 15.5), (16.5, 17.0), (11.5, 16.5), (7.0, 14.0), (3.0, 11.5), (0.0, 7.0),
])
outline = Path(catchment)
x, y = catchment.T
area = 0.5 * abs(x @ np.roll(y, -1) - y @ np.roll(x, -1))   # shoelace formula
town = np.array([4.5, 4.0])

def rain(x, y):
    """Storm total in mm: rising from the south-west toward the ridge, plus one storm cell."""
    ridge = 20 + 35 * (0.6 * (x - 0.5) / 23.5 + 0.4 * y / 17.0)
    cell = 12 * np.exp(-((x - 13) ** 2 + (y - 9) ** 2) / (2 * 3.5 ** 2))
    return ridge + cell

rng = np.random.default_rng(1911)
gauges = []
while len(gauges) < 5:                                    # five by the town, where the observers live
    p = town + 1.2 * rng.standard_normal(2)
    if outline.contains_point(p) and all(np.hypot(*(p - q)) > 0.8 for q in gauges):
        gauges.append(p)
lo, hi = catchment.min(axis=0), catchment.max(axis=0)
while len(gauges) < 12:                                   # seven over the rest, at least 3 km apart
    p = rng.uniform(lo, hi)
    if outline.contains_point(p) and all(np.hypot(*(p - q)) > 3.0 for q in gauges):
        gauges.append(p)
gauges = np.array(gauges)
P = rain(*gauges.T) + rng.normal(0, 1.0, len(gauges))   # readings in mm, 1 mm gauge noise

# the true mean: the field averaged over a fine grid inside the outline
dx = 0.03                                                 # km
X, Y = np.meshgrid(np.arange(lo[0], hi[0], dx) + dx / 2, np.arange(lo[1], hi[1], dx) + dx / 2)
grid = np.column_stack([X.ravel(), Y.ravel()])
inside = outline.contains_points(grid)
true_mean = rain(*grid[inside].T).mean()

def draw_catchment(ax):
    ax.fill(*catchment.T, color=INK, alpha=0.04, lw=0)
    ax.plot(*np.vstack([catchment, catchment[:1]]).T, color=INK, lw=1.2)
    ax.set(xlim=(lo[0] - 0.8, hi[0] + 0.8), ylim=(lo[1] - 0.8, hi[1] + 0.8), aspect="equal",
           xticks=range(0, 25, 5), yticks=range(0, 20, 5), xlabel="x / km", ylabel="y / km")

NUDGE = {0: (-4, 6), 3: (-16, 5), 4: (5, -13)}            # three town labels move aside, in points
def label_readings(ax):
    for k, ((xg, yg), p) in enumerate(zip(gauges, P)):
        ax.annotate(f"{p:.0f}", (xg, yg), xytext=NUDGE.get(k, (5, 4)), textcoords="offset points",
                    color=INK, path_effects=[withStroke(linewidth=3, foreground="white")])

fig, ax = plt.subplots()
draw_catchment(ax)
ax.plot(*gauges.T, "o", color=INK, ms=6)
label_readings(ax)
ax.text(*(town + [0.6, -2.9]), "town", color=MUTED, ha="center")
plt.show()
print(f"catchment: {area:.1f} km², {len(gauges)} gauges, readings {P.min():.0f} to {P.max():.0f} mm")
print(f"plain average of the readings: {P.mean():.1f} mm")

from scipy.spatial import Voronoi

# The exact boundaries come from scipy.spatial.Voronoi; "See it in code" explains the four far points.
ghosts = (lo + hi) / 2 + 240 * np.array([(1, 0), (0, 1), (-1, 0), (0, -1)])
vor_closed = Voronoi(np.vstack([gauges, ghosts]))
ridges = [(i, j, vor_closed.vertices[r]) for (i, j), r in zip(vor_closed.ridge_points, vor_closed.ridge_vertices)
          if i < len(gauges) and j < len(gauges)]

def clip_to_disk(seg, c, r):
    """The part of segment seg within distance r of point c, or None: the boundary burnt so far."""
    a, d = seg[0], seg[1] - seg[0]
    qa, qb, qc = d @ d, 2 * d @ (a - c), (a - c) @ (a - c) - r * r
    disc = qb * qb - 4 * qa * qc
    if disc <= 0:
        return None
    t1, t2 = max((-qb - np.sqrt(disc)) / (2 * qa), 0), min((-qb + np.sqrt(disc)) / (2 * qa), 1)
    return a + np.outer([t1, t2], d) if t1 < t2 else None

# distance to the nearest gauge on a coarse grid over the view: the burnt area and its front
xv, yv = np.meshgrid(np.linspace(lo[0] - 0.8, hi[0] + 0.8, 400), np.linspace(lo[1] - 0.8, hi[1] + 0.8, 300))
d_view = np.min(np.hypot(xv[..., None] - gauges[:, 0], yv[..., None] - gauges[:, 1]), axis=-1)

fig, axes = plt.subplots(1, 3, figsize=(8, 2.6), sharex=True, sharey=True, layout="constrained")
for ax, r in zip(axes, [1.5, 3.0, 5.0]):
    draw_catchment(ax)
    ax.contourf(xv, yv, d_view, levels=[0, r], colors=[MUTED], alpha=0.3)
    ax.contour(xv, yv, d_view, levels=[r], colors=[MUTED], linewidths=1)
    for i, j, seg in ridges:
        part = clip_to_disk(seg, gauges[i], r)
        if part is not None:
            ax.plot(*part.T, color=ACCENT, lw=1.6)
    ax.plot(*gauges.T, "o", color=INK, ms=3)
    ax.text(0.03, 0.97, f"r = {r:.1f} km", transform=ax.transAxes, va="top",
            path_effects=[withStroke(linewidth=3, foreground="white")])
    ax.set_ylabel("y / km" if ax is axes[0] else "")
plt.show()

# distance from every grid point inside the catchment to every gauge, by broadcasting
pts = grid[inside]
dist = np.hypot(pts[:, None, 0] - gauges[:, 0], pts[:, None, 1] - gauges[:, 1])
nearest = dist.argmin(axis=1)
N = np.bincount(nearest, minlength=len(gauges))
A_grid = N * dx * dx                                      # points times the area of one grid square
share = A_grid / A_grid.sum()
print(f"grid: {len(pts):,} points inside the outline, {dx} km apart, {dx * dx * 1e6:.0f} m² each")
for p, s in sorted(zip(P, share), key=lambda ps: -ps[1])[:3]:
    print(f"  gauge reading {p:4.1f} mm: {100 * s:4.1f} % of the area")
print(f"the five town gauges: {100 * share[:5].sum():.0f} % of the area, {100 * 5 / 12:.0f} % of the votes")

from scipy.spatial import Delaunay

tri = Delaunay(gauges)
vor = Voronoi(gauges)                                     # no far points: the diagram of the twelve alone

def circumcircle(a, b, c):
    """Center and radius of the circle through three points."""
    z = np.linalg.solve(2 * np.array([b - a, c - a]), [b @ b - a @ a, c @ c - a @ a])
    return z, np.hypot(*(z - a))

circles = [circumcircle(*gauges[s]) for s in tri.simplices]
centers = np.array([z for z, r in circles])
radii = np.array([r for z, r in circles])

# two circles worth showing: centered in the catchment, not too large, as far apart as possible
ok = np.flatnonzero(outline.contains_points(centers) & (radii < 6))
pairs = [(k, m) for k in ok for m in ok if k < m]
shown = max(pairs, key=lambda km: np.hypot(*(centers[km[0]] - centers[km[1]])))

def segments_cross(a, b, c, d):
    side = lambda p, q, r: np.sign((q - p)[0] * (r - p)[1] - (q - p)[1] * (r - p)[0])
    return side(a, b, c) != side(a, b, d) and side(c, d, a) != side(c, d, b)

fig, ax = plt.subplots()
draw_catchment(ax)
for i, j, seg in ridges:
    ax.plot(*seg.T, color=ACCENT, lw=1.6)
ax.triplot(*gauges.T, tri.simplices, color=SECOND, lw=1.2)
for k in shown:
    ax.add_patch(plt.Circle(centers[k], radii[k], fill=False, color=MUTED, lw=1, ls="--"))
    ax.plot(*centers[k], "o", color=MUTED, ms=4)
ax.plot(*gauges.T, "o", color=INK, ms=6, zorder=3)
plt.show()
crossing = sum(segments_cross(gauges[i], gauges[j], *seg) for i, j, seg in ridges)
print(f"{len(tri.simplices)} Delaunay triangles, {len(vor.vertices)} Voronoi cell corners")
print(f"{crossing} of {len(ridges)} Delaunay edges cross their cell edge")

from scipy.spatial import ConvexHull

# for each triangle: distance from its circle's center to the nearest gauge not on it, over the radius
ratio = min(np.delete(np.hypot(*(gauges - z).T), s).min() / r for (z, r), s in zip(circles, tri.simplices))
n, h = len(gauges), len(ConvexHull(gauges).vertices)
edges = {tuple(sorted((s[a], s[b]))) for s in tri.simplices for a, b in [(0, 1), (1, 2), (2, 0)]}
gap = max(np.hypot(*(centers - v).T).min() for v in vor.vertices)
print(f"empty circle: smallest ratio of nearest other gauge to radius = {ratio:.2f}")
print(f"n = {n}, on the hull h = {h}")
print(f"triangles {len(tri.simplices):3d}   2n - 2 - h = {2 * n - 2 - h}")
print(f"edges     {len(edges):3d}   3n - 3 - h = {3 * n - 3 - h}")
print(f"largest distance from a cell corner to its circumcenter: {gap:.1e} km")

def angles(a, b, c):
    """The three angles of a triangle, in degrees."""
    out = []
    for p, q, r in [(a, b, c), (b, c, a), (c, a, b)]:
        u, v = q - p, r - p
        out.append(np.degrees(np.arccos(u @ v / (np.hypot(*u) * np.hypot(*v)))))
    return out

def cross(o, p, q):
    return (p - o)[0] * (q - o)[1] - (p - o)[1] * (q - o)[0]

flips = []
for k, s in enumerate(tri.simplices):
    for m, nb in enumerate(tri.neighbors[k]):
        if nb <= k:
            continue                                      # each interior edge once
        e1, e2 = np.delete(s, m)                          # the shared edge
        o1, o2 = s[m], (set(tri.simplices[nb]) - {e1, e2}).pop()   # the two gauges facing it
        # the flip is allowed only if the new diagonal o1-o2 crosses the old one
        if np.sign(cross(gauges[o1], gauges[o2], gauges[e1])) == np.sign(cross(gauges[o1], gauges[o2], gauges[e2])):
            continue
        before = min(angles(*gauges[[o1, e1, e2]]) + angles(*gauges[[o2, e1, e2]]))
        after = min(angles(*gauges[[e1, o1, o2]]) + angles(*gauges[[e2, o1, o2]]))
        flips.append((before - after, before, after, o1, o2, e1, e2))
_, before, after, o1, o2, e1, e2 = max(flips)

quad = [o1, e1, o2, e2]
fig, axes = plt.subplots(1, 2, figsize=(7.5, 4.2), sharex=True, sharey=True, layout="constrained")
for ax, pair, color, worst in [(axes[0], [(o1, e1, e2), (o2, e1, e2)], SECOND, before),
                               (axes[1], [(e1, o1, o2), (e2, o1, o2)], MUTED, after)]:
    for t in pair:                                        # Delaunay pair in SECOND, the flipped one in MUTED
        ax.fill(*gauges[list(t)].T, fill=False, ec=color, lw=1.6)
        z, r = circumcircle(*gauges[list(t)])
        ax.add_patch(plt.Circle(z, r, fill=False, color=MUTED, lw=1, ls="--"))
    ax.plot(*gauges[quad].T, "o", color=INK, ms=6, zorder=3)
    for t in pair:                                        # ring the gauge of the four that a circle holds
        z, r = circumcircle(*gauges[list(t)])
        swallowed = [g for g in quad if g not in t and np.hypot(*(gauges[g] - z)) < r]
        ax.plot(*gauges[swallowed].T, "o", mfc="none", mec=INK, mew=1.4, ms=15, zorder=4)
    ax.text(0.03, 0.03, f"smallest angle {worst:.0f}°", transform=ax.transAxes,
            path_effects=[withStroke(linewidth=3, foreground="white")])
    ax.set(aspect="equal", xlabel="x / km")
axes[0].set(xlim=(gauges[quad, 0].min() - 1.5, gauges[quad, 0].max() + 1.5),
            ylim=(gauges[quad, 1].min() - 2.5, gauges[quad, 1].max() + 1.5), ylabel="y / km")
plt.show()
print(f"smallest angle: Delaunay pair {before:.0f}°, flipped pair {after:.0f}°")

from shapely.geometry import Polygon

ghosts = (lo + hi) / 2 + 240 * np.array([(1, 0), (0, 1), (-1, 0), (0, -1)])   # far outside, one per side
vor_closed = Voronoi(np.vstack([gauges, ghosts]))
basin = Polygon(catchment)
cells = []
for k in range(len(gauges)):
    region = vor_closed.regions[vor_closed.point_region[k]]
    assert -1 not in region                               # every real cell is closed
    cells.append(Polygon(vor_closed.vertices[region]).intersection(basin))
A = np.array([c.area for c in cells])
thiessen = (A * P).sum() / A.sum()
tri = Delaunay(gauges)

print(f"plain average      {P.mean():5.1f} mm   ({P.mean() - true_mean:+.1f} mm)")
print(f"Thiessen average   {thiessen:5.1f} mm   ({thiessen - true_mean:+.1f} mm)")
print(f"true mean          {true_mean:5.1f} mm")
print(f"cell areas sum to {A.sum():.2f} km², the catchment is {basin.area:.2f} km²")
print(f"largest difference from the grid count: {100 * np.max(np.abs(A_grid - A) / A):.2f} %")

fig, ax = plt.subplots(figsize=(7.5, 5))
draw_catchment(ax)
norm = plt.Normalize(15, 55)
for c, p in zip(cells, P):
    ax.fill(*c.exterior.xy, color=RAIN(norm(p)), lw=0)
    ax.plot(*c.exterior.xy, color=ACCENT, lw=1.4)
for line in ax.triplot(*gauges.T, tri.simplices, color=SECOND, lw=1.2):   # white halo over the shading
    line.set_path_effects([withStroke(linewidth=3, foreground="white")])
ax.plot(*gauges.T, "o", color=INK, ms=6, zorder=3)
label_readings(ax)
ax.set_ylim(top=hi[1] + 3.2)                              # a band above the map for the three averages
ax.text(0.02, 0.98, f"plain average   {P.mean():.1f} mm\narea-weighted   {thiessen:.1f} mm\n"
        f"true mean       {true_mean:.1f} mm", transform=ax.transAxes, va="top", family="monospace")
fig.colorbar(plt.cm.ScalarMappable(norm, RAIN), ax=ax, label="rainfall / mm", shrink=0.8)
plt.show()