Street networks, routing and isochrones

"""Dijkstra's algorithm on six junctions and seven streets, run twice: once by length in metres, once by driving time
in seconds. The street from C to D is a 60 km/h road; every other street is a 20 km/h lane."""
import heapq

STREETS = [  # (junction, junction, length in m, speed in km/h)
    ("A", "B", 900, 20), ("A", "C", 600, 20), ("B", "E", 500, 20), ("B", "F", 1_100, 20),
    ("C", "D", 1_800, 60), ("D", "F", 700, 20), ("E", "F", 800, 20),
]


def graph(cost):
    """Each junction's neighbours with the cost of the street to them, both ways."""
    g = {}
    for a, b, metres, kmh in STREETS:
        w = metres if cost == "length" else round(metres / (kmh / 3.6))  # seconds at that speed
        g.setdefault(a, []).append((b, w))
        g.setdefault(b, []).append((a, w))
    return g


def dijkstra(g, start):
    """Settle the nearest junction not yet settled, then try every street out of it."""
    best, came_from, settled, queue = {start: 0}, {}, [], [(0, start)]
    while queue:
        d, node = heapq.heappop(queue)
        if node in settled:
            continue  # an older, longer entry for a junction already settled
        settled.append(node)
        for nxt, w in g[node]:
            if d + w < best.get(nxt, float("inf")):
                best[nxt], came_from[nxt] = d + w, node
                heapq.heappush(queue, (d + w, nxt))
    return best, came_from, settled


def path(came_from, end):
    out = [end]
    while out[-1] in came_from:
        out.append(came_from[out[-1]])
    return "-".join(reversed(out))


metres = {(a, b): m for a, b, m, _ in STREETS} | {(b, a): m for a, b, m, _ in STREETS}
seconds = {(a, b): round(m / (v / 3.6)) for a, b, m, v in STREETS} | {(b, a): round(m / (v / 3.6)) for a, b, m, v in STREETS}
for cost, unit in (("length", "m"), ("time", "s")):
    best, came_from, settled = dijkstra(graph(cost), "A")
    print(f"by {cost}: settled in the order " + ", ".join(f"{n} {best[n]:,} {unit}" for n in settled))
    route = path(came_from, "F").split("-")
    legs = list(zip(route, route[1:]))
    print(f"  A to F: {'-'.join(route)}, {sum(metres[s] for s in legs):,} m in {sum(seconds[s] for s in legs)} s")

by length: settled in the order A 0 m, C 600 m, B 900 m, E 1,400 m, F 2,000 m, D 2,400 m
  A to F: A-B-F, 2,000 m in 360 s
by time: settled in the order A 0 s, C 108 s, B 162 s, D 216 s, E 252 s, F 342 s
  A to F: A-C-D-F, 3,100 m in 342 s

"""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 track's town as street networks, in metres of UTM zone 43N. For walking: a junction every 100 m, a river at
1,950 m north of the town's south edge that only three bridges cross, and a park with no streets inside; every street
can be walked both ways. For driving: the same streets as one-way or two-way edges with a speed each."""
import networkx as nx
import shapely
from town import E0, N0

W, H, STEP = 8_000, 6_000, 100
RIVER = N0 + 1_950  # between the rows of junctions at 1,900 and 2,000 m north
BRIDGES = (E0 + 1_000, E0 + 3_500, E0 + 6_000)
PARK = (E0 + 4_400, N0 + 4_200, E0 + 5_400, N0 + 4_800)  # west, south, east, north


def streets():
    """Every street between two neighbouring junctions, as ((x1, y1), (x2, y2))."""
    out = []
    for x in range(E0, E0 + W + 1, STEP):
        for y in range(N0, N0 + H + 1, STEP):
            for end in ((x + STEP, y), (x, y + STEP)):
                if end[0] > E0 + W or end[1] > N0 + H:
                    continue
                if y < RIVER < end[1] and x not in BRIDGES:
                    continue  # no bridge here
                mx, my = (x + end[0]) / 2, (y + end[1]) / 2
                if PARK[0] < mx < PARK[2] and PARK[1] < my < PARK[3]:
                    continue  # inside the park
                out.append(((x, y), end))
    return out


def walk_network():
    g = nx.Graph()
    g.add_edges_from(streets(), length=STEP)
    return g


def nearest_junction(g):
    """A function that gives the junction nearest to a point: the corner of its block, or the park's edge for a
    point inside the park, found with a spatial index of the junctions."""
    nodes = list(g)
    tree = shapely.STRtree(shapely.points(nodes))

    def nearest(x, y):
        node = (round(x / STEP) * STEP, round(y / STEP) * STEP)
        return node if node in g else nodes[tree.nearest(shapely.Point(x, y))]

    return nearest


def speed(a, b):
    """km/h: the ring road around the town's edge 60, Station Road (x = 503,500) 40, every other street 20."""
    if a[0] == b[0] in (E0, E0 + W) or a[1] == b[1] in (N0, N0 + H):
        return 60
    if a[0] == b[0] == E0 + 3_500:
        return 40
    return 20


def drive_network():
    """Directed: the streets of the market blocks (x 503,000 to 504,000, y 2,002,500 to 2,003,500) are one-way,
    eastwards and northwards; every other street runs both ways."""
    g = nx.DiGraph()
    for a, b in streets():
        minutes = STEP / (speed(a, b) * 1000 / 60)
        g.add_edge(a, b, length=STEP, minutes=minutes)
        market = E0 + 3_000 <= min(a[0], b[0]) and max(a[0], b[0]) <= E0 + 4_000 and N0 + 2_500 <= min(a[1], b[1]) and max(a[1], b[1]) <= N0 + 3_500
        if not market:
            g.add_edge(b, a, length=STEP, minutes=minutes)
    return g

"""Where can you walk in 10 minutes from clinic C2? Dijkstra's algorithm over the town's streets finds every junction
within 800 m on foot, at 80 m a minute; the streets between them, widened by half a block on each side, make the
walk's service area. Beside it, the 800 m circle that a buffer around C2 would draw."""
import math

import networkx as nx
import shapely
from streets import STEP, nearest_junction, walk_network
from town import CLINICS, E0, N0, households

SPEED = 80  # metres a minute on foot, 4.8 km/h
REACH = 10 * SPEED

g = walk_network()
junction = nearest_junction(g)
c2 = (E0 + CLINICS["C2"][0], N0 + CLINICS["C2"][1])  # C2 stands on a junction
on_foot = nx.single_source_dijkstra_path_length(g, c2, weight="length")  # metres from C2 to every junction
reached = {n for n, d in on_foot.items() if d <= REACH}
streets = [shapely.LineString([a, b]) for a, b in g.edges() if a in reached and b in reached]
walk = shapely.union_all(shapely.buffer(streets, STEP / 2))
circle = shapely.Point(c2).buffer(REACH, quad_segs=16)
print(f"junctions within a 10-minute walk of C2: {len(reached)}")
shared = walk.intersection(circle).area  # widened by half a block, the walk's four tips reach a little past the circle
print(f"the walk's service area: {walk.area / 1e6:.3f} km², {shared / 1e6:.3f} km² of it inside the 800 m circle")
print(f"the circle: {circle.area / 1e6:.3f} km², of which the walk covers {shared / circle.area:.1%}")

inside = [n for n in g if math.dist(n, c2) <= REACH]
far = max(inside, key=lambda n: (on_foot[n], n))
dx, dy = far[0] - c2[0], far[1] - c2[1]
print(f"the circle holds {len(inside)} junctions; the farthest on foot, {abs(dx)} m {'east' if dx > 0 else 'west'} and {abs(dy)} m "
      f"{'north' if dy > 0 else 'south'} of C2, is {math.dist(far, c2):.0f} m away in a straight line and {on_foot[far]:,} m on foot: {on_foot[far] / SPEED:.1f} minutes")


def walk_metres(x, y):
    """Straight to the nearest junction, then along the streets."""
    node = junction(x, y)
    return math.dist((x, y), node) + on_foot[node]


homes = households()
by_circle = sum(1 for _, x, y in homes if math.dist((x, y), c2) <= REACH)
by_walk = sum(1 for _, x, y in homes if walk_metres(x, y) <= REACH)
print(f"households: {by_circle} within 800 m in a straight line, {by_walk} within a 10-minute walk")

junctions within a 10-minute walk of C2: 144
the walk's service area: 1.421 km², 1.405 km² of it inside the 800 m circle
the circle: 2.007 km², of which the walk covers 70.0%
the circle holds 189 junctions; the farthest on foot, 600 m east and 500 m north of C2, is 781 m away in a straight line and 1,100 m on foot: 13.8 minutes
households: 160 within 800 m in a straight line, 87 within a 10-minute walk