GIS & Geospatial Data (QGIS, GeoPandas, PostGIS) Module 10 – Spatial analysis and statistics
Proximity, service areas and suitability
Measure access to clinics in a straight line, on foot and by car, draw Voronoi catchments, and rank sites for a new clinic with weights you then test.
What you will learn
- Measure access to facilities by straight-line distance, walking distance and driving time, and compare the wards each flags
- Build Voronoi catchments for a set of facilities and state the assumptions they make
- Combine criteria in a weighted overlay with hard constraints and a documented table of weights
- Test how sensitive a suitability ranking is to its weights
Before you start
What you’ll learn · Before you start
What you will learn
- Measure access to facilities by straight-line distance, walking distance and driving time, and compare the wards each flags
- Build Voronoi catchments for a set of facilities and state the assumptions they make
- Combine criteria in a weighted overlay with hard constraints and a documented table of weights
- Test how sensitive a suitability ranking is to its weights
Before you start
On this page
In the track’s town, 3,862 of the 5,000 households live within 1,200 m of a clinic in a straight line, yet only 2,549 can walk to one in 15 minutes, which is the same 1,200 m along the streets. Which households count as served depends on how the distance is measured, and every decision built on the count inherits that choice.
This lesson measures access three ways and compares the wards each measure flags, draws the catchments that a nearest-facility rule implies, and then searches the town for a ninth clinic with a weighted overlay: constraints, factors, a table of weights and a test of how much the answer depends on them.
Access depends on how you measure the distance
Three measures answer “how far is the nearest clinic?” for every household:
- Straight line. The distance as the crow flies, in metres of a projected CRS. It needs only two layers of points, and no real trip is shorter.
- Walking distance. The distance along the streets, from one run of Dijkstra's algorithm: A way to find the shortest routes from a start in a graph with weights of 0 or more: settle the nearest junction not yet settled, then check whether going through it shortens the way to each neighbour. More in “Street networks, routing and isochrones” that starts at all eight clinics at once, as in the routing lesson. A walk goes round blocks and crosses the river only at its three bridges.
- Driving time. Minutes on the drive network: 60 km/h on the ring road, 40 on Station Road, 20 elsewhere, and the market’s streets one-way. The search runs on the network reversed, so that it finds each junction’s drive to a clinic.
Each measure needs a limit. A 15-minute walk at 80 m a minute is 1,200 m, so the straight line gets the same 1,200 m, and the drive gets 4 minutes. A ward where fewer than half of the households are within the limit is flagged as underserved. The page shows the program’s two helpers first: the town, and its streets as graphs in plain Python:
"""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 plain Python: 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. Walking uses every street both
ways; driving uses speeds (60 km/h on the ring road, 40 on Station Road, 20 elsewhere) and the market's one-way
streets. Dijkstra's algorithm from every clinic at once gives each junction its distance or time to the nearest one."""
import heapq
import math
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)), always eastwards or northwards."""
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 speed(a, b):
"""km/h of the street from a to b."""
if a[0] == b[0] in (E0, E0 + W) or a[1] == b[1] in (N0, N0 + H):
return 60 # the ring road round the town's edge
if a[0] == b[0] == E0 + 3_500:
return 40 # Station Road
return 20
def one_way(a, b):
"""The streets of the market blocks run only eastwards and northwards."""
return 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
def walk_graph():
"""{junction: [(neighbour, metres), ...]}, every street both ways."""
g = {}
for a, b in streets():
g.setdefault(a, []).append((b, STEP))
g.setdefault(b, []).append((a, STEP))
return g
def drive_graph_to_clinics():
"""{junction: [(previous junction, minutes), ...]}: the drive network reversed, so that Dijkstra's algorithm from
the clinics gives every junction its drive TO the nearest clinic, one-way streets included."""
g = {}
for a, b in streets():
minutes = STEP / (speed(a, b) * 1_000 / 60)
g.setdefault(b, []).append((a, minutes)) # a drive from a to b, reversed
if not one_way(a, b):
g.setdefault(a, []).append((b, minutes)) # and from b to a
return g
def nearest_sources(g, sources):
"""Dijkstra's algorithm from all sources at once: {junction: (cost to the nearest source, that source)}."""
best = {}
heap = [(0.0, name, node) for name, node in sources.items()]
heapq.heapify(heap)
while heap:
cost, name, node = heapq.heappop(heap)
if node in best:
continue
best[node] = (cost, name)
for nxt, w in g.get(node, ()):
if nxt not in best:
heapq.heappush(heap, (cost + w, name, nxt))
return best
def snapper(g):
"""A function that gives the junction nearest to a point: the corner of its block, or for a point inside the park
the nearest junction that has a street."""
nodes = list(g)
def nearest(x, y):
node = (round(x / STEP) * STEP, round(y / STEP) * STEP)
return node if node in g else min(nodes, key=lambda n: math.dist(n, (x, y)))
return nearest"""Access to the 8 clinics measured three ways for every household: in a straight line, on foot along the streets and
by car on the drive network. Each measure has its limit; a ward where fewer than half the households are within it
is flagged as underserved."""
import math
from streets import drive_graph_to_clinics, nearest_sources, snapper, walk_graph
from town import CLINICS, E0, N0, households, ward_bounds
WALK = 80 # metres a minute on foot
LINE_LIMIT = 1_200 # metres in a straight line: a 15-minute walk as the crow flies
WALK_LIMIT = 15 # minutes on foot along the streets
DRIVE_LIMIT = 4 # minutes by car
clinics = {name: (E0 + x, N0 + y) for name, (x, y) in CLINICS.items()}
walk_net = walk_graph()
on_foot = nearest_sources(walk_net, clinics) # {junction: (metres, clinic)}
by_car = nearest_sources(drive_graph_to_clinics(), clinics) # {junction: (minutes, clinic)}
junction = snapper(walk_net)
wards = ward_bounds()
def ward_of(x, y):
return next(w for w, (x0, y0, x1, y1) in wards.items() if x0 <= x < x1 and y0 <= y < y1)
access = {} # household: (ward, straight-line metres, walking minutes, driving minutes)
for h, x, y in households():
j = junction(x, y)
to_street = math.dist((x, y), j) # walked (or driven at 20 km/h) to the nearest junction
line = min(math.dist((x, y), c) for c in clinics.values())
walk = (to_street + on_foot[j][0]) / WALK
drive = to_street / (20_000 / 60) + by_car[j][0]
access[h] = (ward_of(x, y), line, walk, drive)
print(f"households: {len(access):,}; within reach of a clinic")
print(f" {LINE_LIMIT:,} m in a straight line: {sum(a[1] <= LINE_LIMIT for a in access.values()):,}")
print(f" {WALK_LIMIT} minutes on foot: {sum(a[2] <= WALK_LIMIT for a in access.values()):,}")
print(f" {DRIVE_LIMIT} minutes by car: {sum(a[3] <= DRIVE_LIMIT for a in access.values()):,}")
print()
print("share of each ward's households within reach (* fewer than half: underserved)")
print(" ward homes straight line on foot by car")
flagged = {"straight line": [], "on foot": [], "by car": []}
for w in wards:
rows = [a for a in access.values() if a[0] == w]
shares = [sum(r[i] <= lim for r in rows) / len(rows) for i, lim in ((1, LINE_LIMIT), (2, WALK_LIMIT), (3, DRIVE_LIMIT))]
for name, share in zip(flagged, shares):
if share < 0.5:
flagged[name].append(w)
cells = "".join(f"{s:>11.1%} {'*' if s < 0.5 else ' '}" for s in shares)
print(f" {w} {len(rows):>5,} {cells}")
for name, ws in flagged.items():
print(f"underserved {name}: {', '.join(ws)}")Output
households: 5,000; within reach of a clinic 1,200 m in a straight line: 3,862 15 minutes on foot: 2,549 4 minutes by car: 3,154 share of each ward's households within reach (* fewer than half: underserved) ward homes straight line on foot by car W01 295 90.2% 64.4% 76.9% W02 696 81.3% 34.2% * 64.7% W03 183 71.0% 50.8% 58.5% W04 1,169 73.1% 53.0% 60.8% W05 204 36.8% * 7.4% * 16.2% * W06 416 92.8% 66.8% 72.8% W07 379 93.7% 68.3% 79.9% W08 340 60.3% 33.5% * 45.3% * W09 126 91.3% 61.9% 76.2% W10 151 2.6% * 0.0% * 0.7% * W11 648 81.2% 62.2% 70.7% W12 393 96.4% 66.7% 79.1% underserved straight line: W05, W10 underserved on foot: W02, W05, W08, W10 underserved by car: W05, W08, W10
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 access.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
The counts differ by more than a thousand households: 3,862 within reach in a straight line, 2,549 on foot and 3,154 by car. So do the flagged wards. The straight line flags W05 and W10; the walk adds W02 and W08; the drive flags W05, W08 and W10. W02, with 696 households, shows the gap: 81.3 % of them are within 1,200 m in a straight line, 34.2 % within a 15-minute walk, and 64.7 % within a 4-minute drive, which covers 1,333 m on a 20 km/h street and more on Station Road. Which measure fits depends on the trip and on the data you have. Answer two questions about yours:
Answer the questions above: the table then marks what fits your case and says why the rest does not.
| Measure | Data it needs | People on foot | Cars and ambulances | In the lesson's town |
|---|---|---|---|---|
| Straight line Measures: metres as the crow flies | Good two layers of points in a projected CRS | Poor ignores blocks and the river: 3,862 within 1.2 km, 2,549 on foot | Poor ignores roads, speeds and one-way streets | Mixed flags only W05 and W10 |
| Walking distance Measures: metres along the streets | Mixed a connected street network with lengths | Good follows the streets and the three bridges | Poor no speeds and no one-way streets | Good flags W02, W05, W08 and W10 |
| Driving time Measures: minutes on roads with speeds | Poor streets with speeds and directions | Poor times a car, not a walk | Good uses the fast roads and the one-way market streets | Mixed flags W05, W08 and W10: 3,154 within 4 minutes |
How the answers decide
- “How do the people in question reach the clinic?” On foot: rules out Driving time and adds a caveat to Straight line. A drive time says how long a car takes, not who can walk there; a straight line counts 1,313 more households than the streets allow.
- “How do the people in question reach the clinic?” By car or ambulance: rules out Walking distance and adds a caveat to Straight line. A car is timed by road speeds and one-way streets, which a walking distance ignores; a straight line ignores both.
- “How do the people in question reach the clinic?” Not known yet: a first screen: adds a caveat to Walking distance and Driving time. For a first screen the straight line is quick and cheap; say that it overstates access, and measure on the network before deciding.
- “Which street data do you have?” No street network: rules out Walking distance and Driving time. Without a street network only the straight line can be computed: report it as an upper bound on access.
- “Which street data do you have?” Streets with lengths, no speeds: rules out Driving time. A drive time needs speeds and directions; without them it is a walking distance with a different unit.
Change the limits
Press Edit and set WALK_LIMIT = 20: the walk now flags only W05 and W10, the straight line’s list. Set it back and
make DRIVE_LIMIT = 5: the drive flags only W10. The list of underserved wards is a choice of limit as much as of
measure, so a report names both.
Checkpoint
Check: Three ways to measure access
Voronoi catchments give every place its nearest clinic
A catchment says which households a facility serves. The simplest rule gives every place to the facility nearest to it in a straight line, and the areas that rule draws form a Voronoi diagram (Thiessen polygons): A division of an area into cells, one around each site, in which every place belongs to the site nearest to it in a straight line. Each edge lies halfway between two sites., also called Thiessen polygons after Alfred Thiessen, who used them in 1911 to average rainfall over large areas from scattered gauges (Thiessen, Monthly Weather Review). Each edge lies halfway between two sites, on the perpendicular bisector of the line that joins them, so a cell can be built by cutting the study area along the bisector between its site and every other site. The program does exactly that for the eight clinics, then counts the households of each cell and checks them against the walk:
"""Voronoi (Thiessen) catchments of the 8 clinics inside the town: each cell is the town's rectangle cut by the
perpendicular bisector between its clinic and every other clinic. Then the cells' households, the households a cell
gives to a clinic that is not the nearest one on foot, and how many of those live across the river from it."""
import math
from streets import RIVER, nearest_sources, snapper, walk_graph
from town import CLINICS, E0, N0, households
TOWN = [(0, 0), (8_000, 0), (8_000, 6_000), (0, 6_000)] # metres from the south-west corner
def clip(polygon, a, b):
"""The part of a convex polygon that is nearer to point a than to point b (Sutherland-Hodgman, one edge)."""
nx, ny = b[0] - a[0], b[1] - a[1] # inside: n . p <= c, the half-plane on a's side of the bisector
c = (b[0] ** 2 + b[1] ** 2 - a[0] ** 2 - a[1] ** 2) / 2
side = [nx * x + ny * y - c for x, y in polygon]
out = []
for i, p in enumerate(polygon):
q, sp, sq = polygon[(i + 1) % len(polygon)], side[i], side[(i + 1) % len(polygon)]
if sp <= 0:
out.append(p)
if sp * sq < 0: # the edge crosses the bisector
t = sp / (sp - sq)
out.append((p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1])))
return out
def area(polygon):
"""The shoelace formula, in square metres."""
return abs(sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(polygon, polygon[1:] + polygon[:1]))) / 2
cells = {}
for name, site in CLINICS.items():
cell = TOWN
for other, other_site in CLINICS.items():
if other != name:
cell = clip(cell, site, other_site)
cells[name] = cell
net = walk_graph()
junction = snapper(net)
on_foot = {name: nearest_sources(net, {name: (E0 + x, N0 + y)}) for name, (x, y) in CLINICS.items()}
homes = [(h, x - E0, y - N0) for h, x, y in households()]
print(f"the 8 cells cover {sum(area(c) for c in cells.values()) / 1e6:.1f} km², the whole town")
print("clinic cell km² households of them nearer another clinic on foot of those, across the river")
moved = across = 0
for name, cell in cells.items():
inside = [(h, x, y) for h, x, y in homes if min(CLINICS, key=lambda c: math.dist((x, y), CLINICS[c])) == name]
elsewhere = far_bank = 0
for h, x, y in inside:
j = junction(E0 + x, N0 + y)
walks = {c: on_foot[c][j][0] for c in CLINICS}
if walks[name] > min(walks.values()):
elsewhere += 1
far_bank += (N0 + y > RIVER) != (N0 + CLINICS[name][1] > RIVER) # the river lies between home and clinic
moved += elsewhere
across += far_bank
print(f" {name} {area(cell) / 1e6:>8.2f} {len(inside):>10,} {elsewhere:>5}{far_bank:>38}")
sizes = [sum(1 for _, x, y in homes if min(CLINICS, key=lambda c: math.dist((x, y), CLINICS[c])) == n) for n in CLINICS]
print(f"households per cell: {min(sizes):,} to {max(sizes):,}, {max(sizes) / min(sizes):.1f} times as many")
print(f"households whose cell is not their nearest clinic on foot: {moved}, of them {across} across the river from their cell's clinic")Output
the 8 cells cover 48.0 km², the whole town clinic cell km² households of them nearer another clinic on foot of those, across the river C1 6.15 572 38 0 C2 5.34 489 37 0 C3 7.84 1,439 22 0 C4 7.02 677 50 23 C5 4.64 454 37 29 C6 7.25 214 3 3 C7 5.91 906 18 12 C8 3.86 249 52 44 households per cell: 214 to 1,439, 6.7 times as many households whose cell is not their nearest clinic on foot: 257, of them 111 across the river from their cell's clinic
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 catchments.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
The 8 clinics' Voronoi catchments and the households in each
Text description of the diagram
The town is a rectangle 8 km wide and 6 km tall. Eight straight-edged cells divide it, one around each clinic, and every edge lies halfway between two clinics, on the perpendicular bisector of the line that joins them. Every place inside a cell is nearer to that cell's clinic, in a straight line, than to any other clinic.
The households in each cell: C1 in the north-west, 572; C2 in the north, 489; C3 in the north-east, 1,439, the most, beside the densest cluster of homes; C4 in the west of the middle, 677; C5 in the middle, 454; C6 in the south-west, 214, the fewest; C7 in the south, 906; C8 in the south-east, 249.
A dashed line crosses the town from west to east at 1,950 m north: the river, which only three bridges cross. The cells of C4, C5, C6, C7 and C8 straddle it, so some households in those cells must walk round to a bridge to reach their cell's clinic.
The cells cover the whole town, 48.0 km², with no gap and no overlap, and their households range from 214 in C6’s cell to 1,439 in C3’s, 6.7 times as many. The streets disagree with them, too: 257 households would reach another clinic sooner on foot than the one their cell gives them. For 111 of them the river lies between home and the cell’s clinic, and 52 live in C8’s cell, which the river cuts in two.
QGIS builds the cells with its Voronoi polygons algorithm (native:voronoipolygons), which can widen the extent by
a buffer region (QGIS).
Shapely’s voronoi_polygons takes extend_to, the geometry whose envelope the diagram must cover, and
ordered=True, added in Shapely 2.1.0 and needing GEOS 3.12 or later, which returns the cells in the order of the
input points (Shapely). Intersect
the cells with the study area when its boundary is not a rectangle.
Used as catchments, Voronoi cells make three assumptions, and each is worth checking against the data:
- Straight lines. People go to the clinic nearest as the crow flies; the town’s 257 households show where the streets disagree.
- Equal clinics. Every clinic is as attractive as the next and has room for everyone its cell sends: here 214 to 1,439 households.
- Nearest only. Nobody passes one clinic to reach another for a doctor, a pharmacy or a shorter wait.
Where one of them fails, build catchments from the trips people make, the nearest clinic along the streets or the clinics they are recorded visiting, and keep Voronoi cells for a first look.
Checkpoint
Check: Voronoi catchments and their assumptions
A weighted overlay turns criteria into one score per cell
Where should a ninth clinic go? A suitability analysis answers such a question for every place at once. On one grid,
each constraint (in a suitability analysis): A yes-or-no rule that rules places in or out before any scoring, such as a park or a setback from a river. A constraint carries no weight, so no score can outweigh it. rules places in or out, each factor scores how good an allowed place
is, and a weighted overlay: Combining factor layers on one grid into a score for every cell: each factor is scaled to a common range, such as 0 to 1, multiplied by its weight, and the results are added. adds the factors into one score per cell
(Malczewski, Progress in Planning). It is
map algebra: Computing a new raster from one or more aligned rasters with a formula applied cell by cell, such as a band ratio, a sum or a comparison with a threshold. Kinds that also read a cell's neighbours or a whole zone are moving-window and zonal operations. More in “Map algebra, NDVI and reclassification” with a table of judgements beside it. The layers usually start as vectors and
distances: polygons are burned into the grid with rasterize: To turn vector features into a raster on a chosen grid by writing a value into the pixels each feature covers, by a stated rule such as pixel centre inside or all touched. More in “Zonal statistics and raster-vector conversion”
(rasterio), and a distance layer comes from a tool
such as gdal_proximity, which measures from every pixel to the nearest target pixel
(GDAL). Step through the search for the ninth clinic:
The drawing shows the whole scene; its 7 steps are listed below.
Every step as text
- One grid. Every layer of the search uses the same grid: 1,200 cells of 200 m over the 8 km by 6 km town, so that a cell means the same place in each layer.
- Constraints. A constraint is yes or no. Cells inside the park and cells whose centre lies within 150 m of the river, a made-up rule of this town, cannot hold a clinic at any score.
- Factors. Each factor becomes a score from 0 to 1: waiting households within a 15-minute walk, divided by the largest count; nearness to a fast road; and distance from the nearest clinic, up to 2 km.
- Weights. The weights say how much each factor counts, 0.6 for demand, 0.25 for the road and 0.15 for spread, and they add up to 1, so every score stays between 0 and 1.
- Ranking. The 1,145 cells that pass are sorted by score. The best, 300 m from the ring road in W04, scores 0.865 and has 444 waiting households within a 15-minute walk.
- Testing the weights. Each weight goes 20 % up and 20 % down in turn. The best cell stays first in all six runs, and 9 or 10 of the ten best cells stay among the ten best.
- When a ranking flips. A ranking that flips under such small changes says that the weights, not the data, chose the site. Then the weights go back to the people who set them, with the runs that flipped.
The weights are judgements, so they go into a table with their reasons, and the table goes beside the result:
| Factor | Scaled from 0 to 1 | Weight | Why |
|---|---|---|---|
| Demand: waiting households within a 15-minute walk | count ÷ 469, the largest | 0.6 | the clinic is for people with none in walking reach |
| Road: metres to the ring road or Station Road | 1 − metres ÷ 1,000, and 0 beyond | 0.25 | ambulances and staff arrive by road |
| Spread: metres to the nearest clinic | metres ÷ 2,000, and 1 beyond | 0.15 | a site beside a clinic duplicates it |
A waiting household is one more than a 15-minute walk from every clinic today. The constraints come first: the park, and a setback of 150 m from the river. The program scores every cell that is left and ranks them:
"""The factor and constraint layers of a site search for a ninth clinic, on a grid of 200 m cells over the town.
Constraints rule cells out; factors score every cell from 0 (worst) to 1 (best); weights say how much each counts."""
import heapq
import math
from streets import PARK, RIVER, nearest_sources, snapper, walk_graph
from town import CLINICS, E0, N0, households, ward_bounds
CELL = 200
REACH = 1_200 # metres on foot: a 15-minute walk at 80 m a minute
RIVER_SETBACK = 150 # metres: this made-up town's rule for building beside its river
WEIGHTS = {"demand": 0.6, "road": 0.25, "spread": 0.15} # the table of weights, with its reasons in the lesson
def cells():
"""The centre of every 200 m cell, in UTM metres."""
return [(E0 + x + CELL // 2, N0 + y + CELL // 2) for y in range(0, 6_000, CELL) for x in range(0, 8_000, CELL)]
def constraint(x, y):
"""Why a cell is ruled out, or None."""
if PARK[0] < x < PARK[2] and PARK[1] < y < PARK[3]:
return "park"
if abs(y - RIVER) < RIVER_SETBACK:
return "river setback"
return None
def factors():
"""{cell: {factor: 0..1}} for every cell that no constraint rules out."""
net = walk_graph()
junction = snapper(net)
sites = {name: (E0 + x, N0 + y) for name, (x, y) in CLINICS.items()}
on_foot = nearest_sources(net, sites)
# Households more than a 15-minute walk from every clinic, grouped by the junction they walk to.
waiting = {}
for _, x, y in households():
j = junction(x, y)
snap = math.dist((x, y), j)
if snap + on_foot[j][0] > REACH:
waiting.setdefault(j, []).append(snap)
raw = {}
for x, y in cells():
if constraint(x, y):
continue
start = junction(x, y)
reached = walk_from(net, start, REACH)
demand = sum(1 for j, d in reached.items() for snap in waiting.get(j, ()) if d + snap <= REACH)
road = min(x - E0, E0 + 8_000 - x, y - N0, N0 + 6_000 - y, abs(x - (E0 + 3_500))) # the ring road, Station Road
spread = min(math.dist((x, y), s) for s in sites.values())
raw[(x, y)] = (demand, road, spread)
most = max(r[0] for r in raw.values())
return {c: {"demand": d / most, "road": max(0.0, 1 - r / 1_000), "spread": min(s, 2_000) / 2_000} for c, (d, r, s) in raw.items()}, raw
def walk_from(net, start, limit):
"""Dijkstra's algorithm from one junction, stopping at `limit` metres: {junction: metres}."""
best = {}
heap = [(0, start)]
while heap:
d, node = heapq.heappop(heap)
if node in best or d > limit:
continue
best[node] = d
for nxt, w in net[node]:
if nxt not in best:
heapq.heappush(heap, (d + w, nxt))
return best
def scores(layers, weights):
"""The weighted overlay: the weighted sum of every cell's factors, with the weights scaled to add up to 1."""
total = sum(weights.values())
return {c: sum(weights[k] / total * v for k, v in f.items()) for c, f in layers.items()}
def ward_of(x, y):
return next(w for w, (x0, y0, x1, y1) in ward_bounds().items() if x0 <= x < x1 and y0 <= y < y1)"""A weighted overlay for a ninth clinic: rule cells out with the constraints, score the rest with three factors and
the table of weights, and rank them. Then the walk to a clinic with the best cell built."""
import math
from overlay import REACH, WEIGHTS, cells, constraint, factors, scores, ward_of
from streets import nearest_sources, snapper, walk_graph
from town import CLINICS, E0, N0, households
grid = cells()
ruled_out = {}
for x, y in grid:
why = constraint(x, y)
if why:
ruled_out[why] = ruled_out.get(why, 0) + 1
print(f"cells of 200 m: {len(grid):,}; ruled out: " + ", ".join(f"{n} by the {why}" for why, n in ruled_out.items()))
layers, raw = factors()
score = scores(layers, WEIGHTS)
ranked = sorted(score, key=lambda c: (-score[c], c))
print(f"cells scored: {len(score):,}; weights: " + ", ".join(f"{k} {v}" for k, v in WEIGHTS.items()))
print("rank cell centre (m from SW) ward waiting homes road m clinic m score")
for i, c in enumerate(ranked[:5], 1):
demand, road, spread = raw[c]
print(f" {i} ({c[0] - E0:>5,}, {c[1] - N0:>5,}) {ward_of(*c)} {demand:>13} {road:>6,} {spread:>8,.0f} {score[c]:.3f}")
best = ranked[0]
net = walk_graph()
junction = snapper(net)
for label, extra in (("8 clinics", {}), ("with the best cell", {"new": best})):
sites = {name: (E0 + x, N0 + y) for name, (x, y) in CLINICS.items()} | {k: junction(*v) for k, v in extra.items()}
on_foot = nearest_sources(net, sites)
reached = sum(math.dist((x, y), junction(x, y)) + on_foot[junction(x, y)][0] <= REACH for _, x, y in households())
print(f"households within a 15-minute walk, {label}: {reached:,}")Output
cells of 200 m: 1,200; ruled out: 40 by the river setback, 15 by the park cells scored: 1,145; weights: demand 0.6, road 0.25, spread 0.15 rank cell centre (m from SW) ward waiting homes road m clinic m score 1 (7,700, 5,300) W04 444 300 1,628 0.865 2 (7,900, 5,100) W04 383 100 1,664 0.840 3 (7,700, 5,100) W04 426 300 1,500 0.832 4 (7,500, 5,300) W04 460 500 1,487 0.825 5 (7,500, 5,500) W04 446 500 1,640 0.819 households within a 15-minute walk, 8 clinics: 2,549 households within a 15-minute walk, with the best cell: 2,993
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 suitability.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
Of 1,200 cells, 55 are ruled out, 40 by the river setback and 15 by the park. The best of the other 1,145 lies 300 m from the ring road in W04 and scores 0.865, with 444 waiting households within a 15-minute walk of it. A clinic there would raise the households within a 15-minute walk of a clinic from 2,549 to 2,993. All five best cells lie in W04, near the ring road in the town’s north-east corner.
Distances that come from a rule
A setback from a river, a road or a protected site is a constraint only when a rule sets it. Take the distance from the regulation that applies to the site, cite it in the table, and never invent one: the 150 m here belongs to this lesson’s made-up town.
Checkpoint
Check: Constraints, factors and weights
A ranking that flips under small changes is not a result
Every weight in the table is a judgement that someone else could make differently, and testing how much a suitability result depends on its weights is a method of its own (Chen, Yu and Khan, Environmental Modelling & Software). The simplest sensitivity analysis: Running an analysis again with its inputs or weights changed, one at a time or together, to see how much the result depends on them. goes one weight at a time: raise it by 20 %, keep the others, scale all three to add up to 1 again, rank the cells and compare; then lower it by 20 %, and do the same for each weight. The program adds a run with equal weights, a different view of what matters, and then searches for the smallest change of each weight that moves the first place:
"""How sensitive is the ranking to its weights? Raise and lower each weight by 20 % (the others keep their values,
and all are scaled to add up to 1 again), then count how many of the ten best cells stay among the ten best. Then
find, in steps of 5 %, the smallest change of each weight that moves the first place, and that takes it out of its
ward."""
from overlay import WEIGHTS, factors, scores, ward_of
from town import E0, N0
CHANGE = 0.2 # raise and lower each weight by this share
layers, _ = factors()
def top(weights, n=10):
s = scores(layers, weights)
return sorted(s, key=lambda c: (-s[c], c))[:n]
def name(cell):
return f"({cell[0] - E0:,}, {cell[1] - N0:,}) in {ward_of(*cell)}"
base = top(WEIGHTS)
wards = sorted({ward_of(*c) for c in base})
print(f"the ten best cells with the table's weights lie in {', '.join(wards)}; the best is {name(base[0])}")
print("change weights (scaled) best cell of the ten best, still there")
variants = [(f"{k} {sign}{CHANGE:.0%}", {**WEIGHTS, k: WEIGHTS[k] * (1 + f)}) for k in WEIGHTS for sign, f in (("+", CHANGE), ("-", -CHANGE))]
variants.append(("equal weights", {k: 1 for k in WEIGHTS}))
for label, w in variants:
total = sum(w.values())
shown = ", ".join(f"{w[k] / total:.2f}" for k in WEIGHTS)
ten = top(w)
print(f" {label:<14} {shown:<19} {name(ten[0]):<23} {len(set(ten) & set(base)):>2} of 10")
print("the smallest change of each weight, in steps of 5 % from -95 % to +900 %, that")
home = ward_of(*base[0])
print(f"weight moves the first place takes it out of {home}")
for k in WEIGHTS:
moves, leaves = {}, {}
for direction, steps in (("up", range(1, 181)), ("down", range(-1, -20, -1))):
for i in steps:
winner = top({**WEIGHTS, k: WEIGHTS[k] * (1 + 0.05 * i)}, 1)[0]
if winner != base[0] and direction not in moves:
moves[direction] = f"{0.05 * i:+.0%}"
if ward_of(*winner) != home:
leaves[direction] = f"{0.05 * i:+.0%}, to {name(winner)}"
break
moved = " or ".join(moves.values()) or "never"
print(f" {k:<7} {moved:<31} {' or '.join(leaves.values()) or 'never'}")Output
the ten best cells with the table's weights lie in W04; the best is (7,700, 5,300) in W04 change weights (scaled) best cell of the ten best, still there demand +20% 0.64, 0.22, 0.13 (7,700, 5,300) in W04 9 of 10 demand -20% 0.55, 0.28, 0.17 (7,700, 5,300) in W04 10 of 10 road +20% 0.57, 0.29, 0.14 (7,700, 5,300) in W04 10 of 10 road -20% 0.63, 0.21, 0.16 (7,700, 5,300) in W04 9 of 10 spread +20% 0.58, 0.24, 0.17 (7,700, 5,300) in W04 10 of 10 spread -20% 0.62, 0.26, 0.12 (7,700, 5,300) in W04 9 of 10 equal weights 0.33, 0.33, 0.33 (7,900, 5,100) in W04 6 of 10 the smallest change of each weight, in steps of 5 % from -95 % to +900 %, that weight moves the first place takes it out of W04 demand +200% or -35% +870%, to (2,700, 4,100) in W02 or -85%, to (3,500, 900) in W10 road +55% or -85% +875%, to (3,500, 4,500) in W02 spread +745% never
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 sensitivity.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
The best cell, (7,700, 5,300), stays first in all six runs, and 9 or 10 of the ten best cells stay among the ten best. Small disagreements about the weights would not move the clinic. Equal weights would: they put (7,900, 5,100) first, 283 m away, and keep 6 of the ten best.
Make the changes larger
Press Edit and set CHANGE = 0.5: halving demand’s weight now puts (7,900, 5,100) first, while the other five runs
keep (7,700, 5,300). The larger the change a ranking survives, the less its weights decide it.
The last table shows how far the answer stands. Demand’s weight must fall by 35 % or road’s rise by 55 % before the first place moves at all, and then only to a neighbouring cell in W04. In this search, a site in another ward comes first only when a weight grows almost tenfold, or when demand loses 85 % of its weight.
Key idea
Publish the constraints with their sources, the weights with their reasons and the sensitivity runs beside the ranking: they belong to the result, because they decide it.
Checkpoint
Check: Testing the weights
Practise: rank sites for a ninth clinic
The exercise is the overlay in small: five candidate sites, three factors and two constraints. Rank the sites with min-max scaling and the table’s weights, then count how many one-at-a-time changes of 20 % keep the winner. The sample tests include a site exactly at the river setback, a factor with one value at every site, and two sites close enough for the weights to decide between them:
Exercise · Medium · Python
Rank candidate sites for a clinic, then test the weights
Write two functions in sites.py. Each candidate site is a dictionary such as {"id": "S1", "waiting": 444, "road_m": 300, "clinic_m": 1628, "river_m": 900, "park": False}: the households waiting for a clinic within a walk of it, the metres to the nearest fast road, the metres to the nearest existing clinic, the metres to the river, and whether it lies in the park.
rank_sites(sites, weights) returns the ids of the sites that pass both constraints, best first:
- A site in the park, or less than 150 m from the river, is ruled out. A site exactly 150 m away passes.
- Scale each factor to 0..1 over the sites that pass, with min-max scaling: more
waitingand moreclinic_mare better (the largest value scores 1), lessroad_mis better (the smallest value scores 1). When every passing site has the same value of a factor, each scores 1 on it. weightsgives a weight for"waiting","road_m"and"clinic_m". Divide each by their sum, so that{"waiting": 6, "road_m": 2.5, "clinic_m": 1.5}ranks exactly as{"waiting": 0.6, "road_m": 0.25, "clinic_m": 0.15}.- The score is the weighted sum of the scaled factors. Higher scores rank first; equal scores go in order of id.
winner_holds(sites, weights, change=0.2) tests the ranking: for each of the three weights, raise it by change (multiply it by 1 + change), then lower it (multiply by 1 - change), keeping the other two as they are. Of these six variants, it returns how many give the same best site as weights itself.
Starter code · sites.py
def rank_sites(sites, weights):
"""The ids of the sites that pass the constraints, best first (see the prompt)."""
return []
def winner_holds(sites, weights, change=0.2):
"""How many of the six one-at-a-time weight changes keep the same best site."""
return 0The sample tests · test_sites.py
from sites import rank_sites, winner_holds
WEIGHTS = {"waiting": 0.6, "road_m": 0.25, "clinic_m": 0.15}
SITES = [
{"id": "S1", "waiting": 444, "road_m": 300, "clinic_m": 1628, "river_m": 900, "park": False},
{"id": "S2", "waiting": 383, "road_m": 100, "clinic_m": 1664, "river_m": 1200, "park": False},
{"id": "S3", "waiting": 500, "road_m": 800, "clinic_m": 900, "river_m": 140, "park": False},
{"id": "S4", "waiting": 120, "road_m": 50, "clinic_m": 2000, "river_m": 600, "park": True},
{"id": "S5", "waiting": 300, "road_m": 400, "clinic_m": 2400, "river_m": 150, "park": False},
]
CLOSE = [
{"id": "A", "waiting": 400, "road_m": 500, "clinic_m": 1500, "river_m": 500, "park": False},
{"id": "B", "waiting": 360, "road_m": 300, "clinic_m": 1500, "river_m": 500, "park": False},
{"id": "C", "waiting": 100, "road_m": 900, "clinic_m": 1500, "river_m": 500, "park": False},
]
def test_constraints_rule_sites_out():
"""leaves out the site in the park and the one 140 m from the river, and keeps the one at exactly 150 m"""
ranked = rank_sites(SITES, WEIGHTS)
assert "S3" not in ranked and "S4" not in ranked
assert "S5" in ranked
def test_the_weighted_ranking():
"""ranks the passing sites by their weighted, min-max scaled factors"""
assert rank_sites(SITES, WEIGHTS) == ["S1", "S2", "S5"]
def test_less_road_distance_is_better():
"""scores the road factor so that the nearest site to a fast road gets 1"""
assert rank_sites(SITES, {"waiting": 0.2, "road_m": 0.6, "clinic_m": 0.2}) == ["S2", "S1", "S5"]
def test_more_distance_from_a_clinic_is_better():
"""scores the spread factor so that the site farthest from a clinic gets 1"""
assert rank_sites(SITES, {"waiting": 0.1, "road_m": 0.1, "clinic_m": 0.8}) == ["S5", "S2", "S1"]
def test_more_waiting_households_is_better():
"""with demand alone, ranks the sites by their waiting households, most first"""
assert rank_sites(SITES, {"waiting": 1, "road_m": 0, "clinic_m": 0}) == ["S1", "S2", "S5"]
def test_weights_are_scaled_to_add_up_to_one():
"""gives the same ranking for weights that differ only by a common factor"""
assert rank_sites(SITES, {"waiting": 6, "road_m": 2.5, "clinic_m": 1.5}) == ["S1", "S2", "S5"]
def test_a_factor_with_one_value():
"""does not divide by zero when every passing site has the same value of a factor"""
assert rank_sites(CLOSE, WEIGHTS) == ["B", "A", "C"]
def test_a_clear_winner_holds():
"""keeps the best site in all six one-at-a-time changes of 20 %"""
assert winner_holds(SITES, WEIGHTS) == 6
def test_a_close_winner_does_not():
"""finds the two changes of 20 % that hand the first place to another site"""
assert winner_holds(CLOSE, WEIGHTS) == 4A hint
Filter first, then scale: the minimum and maximum of each factor come from the sites that pass the constraints only. For a factor where less is better, the scaled value is (largest - value) / (largest - smallest). In winner_holds, build each variant with dict(weights) and change one key, then compare rank_sites(sites, variant)[0] with the best site of the original weights.
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
GATE Geomatics Engineering
The GATE 2027 syllabus for Geomatics Engineering (GE) lists buffer, overlay, interpolation and network analysis as basic spatial analysis in its common section (GATE 2027 GE syllabus); this module works through each of them on one town.
Then test yourself on the four ideas, from the three measures of access to the test of the weights. The session starts where your answers are weakest:
Practice that adapts to you
Practise this lesson
12 questions about 4 ideas of this lesson. Each answer is explained at once, and what you get right is kept in this browser.
Interview questions
Warm-up (fresher): What is the difference between a constraint and a factor in a suitability analysis?
A constraint decides whether a place may be used at all. Inside a park, or within a setback that a rule sets from a river, a site is ruled out whatever else it offers: the constraint is a yes-or-no mask applied before any scoring. A factor measures how good an allowed place is, such as the households it would bring within a walk or its distance from a fast road; factors are scaled to a common range, usually 0 to 1, weighted and added up. Never give a constraint a weight, because a high enough score on the other factors could then outweigh a rule that must not be broken. In the lesson’s search, two constraints removed 55 of 1,200 cells before three weighted factors scored the other 1,145.
Two longer questions put catchments and weights in front of a district and a council. The module’s bank gives what each tests, a first step and the follow-ups:
Model answers open with any pass: ₹79 for 1 week ($2.99), one payment that never renews.
-
Fresher to mid
A district plans to staff its 8 clinics in proportion to the households in each clinic's Voronoi cell, which range from 214 to 1,439. Would you accept the plan as it stands? What would you check, and which catchments would you build instead?
-
Mid to senior
Your weighted overlay puts a new clinic in the north-east ward. A councillor says road access deserves far more weight and asks you to raise it until a site in their ward comes first. How do you respond, and what do you show the council?
Flashcards · 8 cards
Recall it, then check
Say the answer to yourself before you look: recalling it is what makes it stick.
-
Which three measures answer "how far is the nearest clinic?", and how far apart were they in the town?
Straight-line distance, walking distance along the streets, and driving time on the drive network. With a limit of 1,200 m or 15 minutes on foot, 3,862 households were in reach in a straight line, 2,549 on foot and 3,154 within a 4-minute drive.
-
How can a single search give each household the walking distance to its nearest clinic?
Run Dijkstra's algorithm once, starting from all eight clinics at the same time. For drive times to a clinic, run it on the network reversed, so that it finds each junction's drive to a clinic, one-way streets included.
-
Why must an access report name both its measure and its limit?
Both change the list of underserved wards: the 15-minute walk flagged W02, W05, W08 and W10, while a 20-minute walk flagged only W05 and W10, the same list as the straight line.
-
What is a Voronoi diagram, and how is one cell built?
The areas that give every place to its nearest site in a straight line, also called Thiessen polygons. Each edge lies on the perpendicular bisector between two sites, so a cell is the study area cut along the bisector between its site and every other site.
-
Which three assumptions does a Voronoi catchment make?
That people travel in straight lines, that every facility is as attractive and has room for everyone its cell sends, and that nobody passes one facility for another. The town's cells held 214 to 1,439 households, and 257 households would reach another clinic sooner on foot.
-
How does a constraint differ from a factor in a suitability analysis?
A constraint is a yes-or-no mask applied before any scoring, such as a park or a setback a rule sets from a river. A factor scores how good an allowed place is, scaled from 0 to 1, weighted and added. A constraint never gets a weight, or a high score elsewhere could outweigh a rule.
-
Where did the weighted overlay put the town's ninth clinic, and what would it change?
Two constraints removed 55 of 1,200 cells. With weights of 0.6 for demand, 0.25 for road and 0.15 for spread, the best remaining cell, in W04 and 300 m from the ring road, scored 0.865; a clinic there would raise the households within a 15-minute walk from 2,549 to 2,993.
-
How does a one-at-a-time sensitivity analysis test a suitability ranking?
Raise one weight by 20 %, rescale all of them to add up to 1, rank again, then lower it by 20 %, and repeat for each weight. The town's best cell stayed first in all six runs; equal weights moved it, 283 m away.
Key takeaways
- Access depends on the measure and its limit: 3,862 households were within 1,200 m in a straight line, 2,549 within a 15-minute walk and 3,154 within a 4-minute drive, and each measure flagged a different list of wards.
- A Voronoi cell holds the places nearest to its facility in a straight line; as a catchment it assumes straight-line trips, equal facilities and nobody passing one facility for another.
- A suitability analysis rules places out with constraints, scales each factor from 0 to 1 and adds the factors with weights recorded in a table with their reasons.
- A one-at-a-time test shows how much a ranking depends on its weights: the town’s best cell held at 20 % either way, and equal weights moved it.
Terms in this lesson
- constraint (in a suitability analysis)
- A yes-or-no rule that rules places in or out before any scoring, such as a park or a setback from a river. A constraint carries no weight, so no score can outweigh it.
- sensitivity analysis
- Running an analysis again with its inputs or weights changed, one at a time or together, to see how much the result depends on them.
- Voronoi diagram (Thiessen polygons)
- A division of an area into cells, one around each site, in which every place belongs to the site nearest to it in a straight line. Each edge lies halfway between two sites.
- weighted overlay
- Combining factor layers on one grid into a score for every cell: each factor is scaled to a common range, such as 0 to 1, multiplied by its weight, and the results are added.
References
- QGIS processing algorithms: vector geometry (Voronoi polygons) (QGIS Project)
- shapely.voronoi_polygons (Shapely developers)
- Precipitation averages for large areas (Thiessen) (Monthly Weather Review (American Meteorological Society))
- GIS-based land-use suitability analysis: a critical overview (Malczewski) (Progress in Planning (Elsevier))
- Spatial sensitivity analysis of multi-criteria weights in GIS-based land suitability evaluation (Chen, Yu and Khan) (Environmental Modelling & Software (Elsevier))
- gdal_proximity (Open Source Geospatial Foundation (GDAL))
- Vector features (rasterio documentation) (rasterio developers)
- GATE 2027 syllabus: Geomatics Engineering (GE) (IIT Madras, organising institute of GATE 2027)
Related tools
Report a problem with this lesson
Kept only in this browser. Your Learn progress