Overlay, dissolve, clip and areal weighting

"""The track's synthetic town in metres of UTM zone 43N (EPSG:32643): 12 square wards of 2 km, 5,000 households
from a seeded generator and 8 clinics. Every lesson of the track asks its questions about this town."""
import random

E0, N0 = 500_000, 2_000_000  # the town's south-west corner: easting and northing in metres
WARD = 2_000  # each ward is a 2 km square
COLS, ROWS = 4, 3  # 4 wards from west to east, 3 from south to north

# Metres east and north of the south-west corner. C2 stands exactly on the edge of W02 and W03, on purpose.
CLINICS = {
    "C1": (1_500, 5_200), "C2": (4_000, 5_000), "C3": (6_500, 4_200), "C4": (2_700, 3_100),
    "C5": (5_200, 2_800), "C6": (900, 1_100), "C7": (5_600, 1_200), "C8": (7_400, 900),
}


def ward_bounds():
    """W01 to W12, north-west first, each as (min x, min y, max x, max y)."""
    out = {}
    for row in range(ROWS):
        for col in range(COLS):
            x0, y0 = E0 + col * WARD, N0 + (ROWS - 1 - row) * WARD
            out[f"W{row * COLS + col + 1:02d}"] = (x0, y0, x0 + WARD, y0 + WARD)
    return out


def households(seed=42, n=5_000):
    """(id, x, y) of n households, most of them around three neighbourhood centres."""
    rng = random.Random(seed)
    centres = [(2_600, 4_300, 700), (5_600, 1_500, 900), (6_900, 4_700, 500)]
    homes = []
    while len(homes) < n:
        if rng.random() < 0.7:
            cx, cy, spread = rng.choice(centres)
            # The sum of three random numbers clusters around its middle, like a bell curve.
            x = cx + spread * 2 * (rng.random() + rng.random() + rng.random() - 1.5)
            y = cy + spread * 2 * (rng.random() + rng.random() + rng.random() - 1.5)
        else:
            x, y = rng.uniform(0, COLS * WARD), rng.uniform(0, ROWS * WARD)
        if 0 < x < COLS * WARD and 0 < y < ROWS * WARD:
            homes.append((f"H{len(homes) + 1:04d}", round(E0 + x, 1), round(N0 + y, 1)))
    return homes


def clinics():
    """(id, x, y) of the 8 clinics."""
    return [(k, E0 + x, N0 + y) for k, (x, y) in CLINICS.items()]

"""The town's layers in EPSG:32643 for overlays: wards, households, clinics and their 1 km buffers, 24 census blocks
(each ward split in two at a seeded line), the clinics' catchments (the area nearer to each clinic than to any other),
and two planning zones whose shared edge was digitised a few millimetres off the ward boundary."""
import random

import geopandas as gpd
import shapely
from town import COLS, E0, N0, ROWS, WARD, clinics, households, ward_bounds

CRS = "EPSG:32643"


def wards():
    bounds = ward_bounds()
    return gpd.GeoDataFrame({"ward": list(bounds)}, geometry=[shapely.box(*b) for b in bounds.values()], crs=CRS)


def homes():
    rows = households()
    return gpd.GeoDataFrame({"household": [h for h, _, _ in rows]}, geometry=shapely.points([(x, y) for _, x, y in rows]), crs=CRS)


def clinic_points():
    rows = clinics()
    return gpd.GeoDataFrame({"clinic": [c for c, _, _ in rows]}, geometry=shapely.points([(x, y) for _, x, y in rows]), crs=CRS)


def rings(distance=1_000):
    c = clinic_points()
    return c.set_geometry(c.buffer(distance, quad_segs=16))


def blocks(seed=7):
    """24 census blocks: each ward cut in two by a north-south line at 500 to 1,500 m from its west edge."""
    rng = random.Random(seed)
    rows = []
    for ward, (x0, y0, x1, y1) in ward_bounds().items():
        cut = x0 + rng.randrange(500, 1_501, 100)
        rows += [(ward + "a", ward, shapely.box(x0, y0, cut, y1)), (ward + "b", ward, shapely.box(cut, y0, x1, y1))]
    b = gpd.GeoDataFrame({"block": [r[0] for r in rows], "ward": [r[1] for r in rows]}, geometry=[r[2] for r in rows], crs=CRS)
    inside = gpd.sjoin(homes(), b, predicate="within").groupby("block").size()
    b["households"] = b["block"].map(inside).fillna(0).astype(int)
    return b


def catchments():
    """The part of the town nearer to each clinic than to any other (Voronoi cells, cut to the town)."""
    c = clinic_points()
    town = shapely.box(E0, N0, E0 + COLS * WARD, N0 + ROWS * WARD)
    cells = shapely.voronoi_polygons(shapely.MultiPoint(list(c.geometry)), extend_to=town)
    rows = []
    for cell in cells.geoms:
        name = c.loc[c.intersects(cell), "clinic"].iloc[0]
        rows.append((name, cell.intersection(town)))
    out = gpd.GeoDataFrame({"clinic": [r[0] for r in rows]}, geometry=[r[1] for r in rows], crs=CRS)
    return out.sort_values("clinic").reset_index(drop=True)


def zones():
    """Two planning zones split along y = 2,004,000, but digitised with the edge a few millimetres off."""
    wobble = [0.003, -0.002, 0.004, -0.001, 0.002]
    edge = [(E0 + i * WARD, N0 + 2 * WARD + d) for i, d in enumerate(wobble)]
    north = shapely.Polygon([(E0, N0 + ROWS * WARD), *edge, (E0 + COLS * WARD, N0 + ROWS * WARD)])
    south = shapely.Polygon([(E0, N0), (E0 + COLS * WARD, N0), *edge[::-1]])
    return gpd.GeoDataFrame({"zone": ["North", "South"]}, geometry=[north, south], crs=CRS)

"""What share of each ward lies within 1 km of a clinic? Buffer the clinics, dissolve the buffers into one shape,
overlay it with the wards, and divide each piece's area by its ward's; compare W12 with the rings overlaid without the
dissolve; then clip the households to the same shape."""
import geopandas as gpd
from layers import homes, rings, wards

w = wards()
cover = rings(1_000).dissolve()  # 8 buffers become one MultiPolygon, overlaps counted once
print(f"dissolved buffers: {len(cover)} row, a {cover.geom_type.iloc[0]} of {cover.area.iloc[0] / 1e6:.3f} km²")

# Each ward cut to the covered area. keep_geom_type=True keeps polygons only: where a buffer meets a ward
# at a single corner, the intersection there is a point, and an area table has no use for it.
pieces = w.overlay(cover[["geometry"]], how="intersection", keep_geom_type=True)
pieces["share"] = pieces.area / (2_000 * 2_000)
share = pieces.set_index("ward")["share"].reindex(w["ward"], fill_value=0)
print(f"{len(pieces)} pieces; share of each ward within 1 km of a clinic:")
for (ward, s) in share.items():
    print(f"  {ward} {s:6.1%}  {'#' * round(s * 20)}")
print(f"the town: {pieces.area.sum() / w.area.sum():.1%} of {w.area.sum() / 1e6:.0f} km²")

# The same overlay with the 8 rings as they are: ground where two rings overlap becomes one piece per ring.
loose = w.overlay(rings(1_000), how="intersection", keep_geom_type=True)
w12 = loose.loc[loose["ward"] == "W12"].area.sum() / (2_000 * 2_000)
print(f"W12 from the 8 rings without the dissolve: {w12:.1%}")

reached = gpd.clip(homes(), cover)
print(f"\nclip(households, cover): {len(reached):,} of 5,000 households, with their own columns {reached.columns.tolist()}")
c2_ring = rings(1_000).query("clinic == 'C2'")
print(f"clip(households, C2's ring alone): {len(gpd.clip(homes(), c2_ring))} households")

dissolved buffers: 1 row, a MultiPolygon of 24.695 km²
11 pieces; share of each ward within 1 km of a clinic:
  W01  59.1%  ############
  W02  55.9%  ###########
  W03  49.3%  ##########
  W04  39.0%  ########
  W05   8.8%  ##
  W06  69.6%  ##############
  W07  75.1%  ###############
  W08  28.4%  ######
  W09  75.5%  ###############
  W10   0.0%  
  W11  55.6%  ###########
  W12  82.9%  #################
the town: 49.9% of 48 km²
W12 from the 8 rings without the dissolve: 85.3%

clip(households, cover): 2,971 of 5,000 households, with their own columns ['household', 'geometry']
clip(households, C2's ring alone): 248 households