5.5. Accessibility and Service Areas#

5.5.1. Learning objectives#

By the end of this section, you should be able to:

  • explain why a network service area differs from a Euclidean buffer;

  • choose the mode, speed, threshold, and cost field before calculating accessibility;

  • compute the nodes reachable within a time or distance cutoff;

  • represent a service area as a geometry and say what each representation overstates;

  • join destinations and population to a service area without double-counting; and

  • test how the answer changes with the threshold and with the data.

Section goal

Turn a facility and a time threshold into a service area that can be defended. A small synthetic street grid with a lake in the middle shows what a buffer, a convex hull, and a network service area each claim, and the same code then runs on Columbia.

5.5.2. Network reachability is not a circle#

A Euclidean buffer answers “what is within this straight-line radius?” A network service area answers “what can be reached through this mode and cost within a threshold?” Barriers, one-way streets, bridges, and disconnected components make the two shapes different.

Choose the mode, origin, threshold, and cost field before calculating accessibility [Geurs and Van Wee, 2004]. For example, a 15-minute walk shed should use a pedestrian graph and a travel-time threshold, not a driving graph or a degree-based buffer.

The speed deserves the same care as the mode. ox.add_edge_speeds [Boeing, 2025] fills each edge from its maxspeed tag, or from the average of the tagged edges, and on a walking network those tagged speeds are the speeds of cars. A travel time built from them would describe a 15-minute drive. The toy walking network below is a 9 by 9 grid of 100 m blocks with a lake removed from the middle, and its two main streets carry a 40 km/h tag. The facility is at the lake’s west shore.

Hide code cell source

import math

import networkx as nx
import osmnx as ox
from shapely.geometry import Point, box

WALK_SPEED_KPH = 4.5


def grid(n=9, spacing=100, skip=lambda i, j: False, artery=lambda i, j, di, dj: False):
    # A grid of two-way streets in a projected CRS, in metres.
    G = nx.MultiDiGraph(crs="EPSG:32615")
    ids = {}
    for i in range(n):
        for j in range(n):
            if not skip(i, j):
                ids[(i, j)] = i * n + j
                G.add_node(ids[(i, j)], x=i * spacing, y=j * spacing)
    for (i, j), a in ids.items():
        for di, dj in ((1, 0), (0, 1)):
            b = ids.get((i + di, j + dj))
            if b is not None:
                tagged = artery(i, j, di, dj)
                attributes = dict(length=float(spacing), highway="primary" if tagged else "residential", oneway=False)
                if tagged:
                    attributes["maxspeed"] = "40"
                for u, v in ((a, b), (b, a)):
                    G.add_edge(u, v, **attributes)
    nx.set_node_attributes(G, ox.stats.count_streets_per_node(G), "street_count")
    return G


lake = lambda i, j: 5 <= i <= 6 and 3 <= j <= 5                       # no streets through the lake
main_streets = lambda i, j, di, dj: (i == 0 and dj == 1) or (j == 0 and di == 1)
streets = grid(9, skip=lake, artery=main_streets)
facility = 4 * 9 + 4                                                    # node (4, 4), at the lake's west shore

imputed = ox.add_edge_travel_times(ox.add_edge_speeds(streets.copy()))
walking = streets.copy()
nx.set_edge_attributes(walking, WALK_SPEED_KPH, "speed_kph")
walking = ox.add_edge_travel_times(walking)

for label, graph in [("speeds imputed by OSMnx", imputed), ("fixed walking speed", walking)]:
    speed = sorted({round(data["speed_kph"], 1) for *_, data in graph.edges(data=True)})
    reached = nx.single_source_dijkstra_path_length(graph, facility, cutoff=5 * 60, weight="travel_time")
    print(f"{label:24} speed {speed} km/h -> {len(reached)} of {graph.number_of_nodes()} nodes within 5 minutes")
speeds imputed by OSMnx  speed [40.0] km/h -> 75 of 75 nodes within 5 minutes
fixed walking speed      speed [4.5] km/h -> 18 of 75 nodes within 5 minutes

With imputed speeds every street is a 40 km/h street, so all 75 nodes are within a five-minute “walk”, while at a documented walking speed only 18 are. The cost field is the assumption that the whole analysis rests on, so set the walking speed explicitly and record it.

Next, compare the network service area with the Euclidean circle of the same time. At 4.5 km/h, five minutes is 375 m of straight-line walking. For comparison, the comfortable walking speeds reported for healthy adults aged 20 to 79 range from about 4.6 to 5.3 km/h [Bohannon, 1997], so 4.5 km/h is a slightly cautious value.

Hide code cell source

nodes, edges = ox.graph_to_gdfs(walking)
reachable = nx.single_source_dijkstra_path_length(walking, facility, cutoff=5 * 60, weight="travel_time")

circle = Point(400, 400).buffer(375)                       # the Euclidean 5-minute walk
inside_circle = set(nodes.index[nodes.geometry.within(circle)])

print(f"nodes inside the 375 m circle:        {len(inside_circle)}")
print(f"nodes reachable in 5 minutes on foot: {len(reachable)}")
print(f"in the circle but not reachable:      {len(inside_circle - set(reachable))}")
print(f"reachable but outside the circle:     {len(set(reachable) - inside_circle)}")
nodes inside the 375 m circle:        39
nodes reachable in 5 minutes on foot: 18
in the circle but not reachable:      21
reachable but outside the circle:     0

The circle contains 39 nodes, but only 18 can be reached, and every reachable node is inside the circle. A network can only shrink a Euclidean area, never enlarge it, because a street route is never shorter than the straight line. The 21 missing nodes are the ones that lie across the lake or beyond the detour around it. The circle would count those people as served.

5.5.3. Reach nodes by distance or time#

nx.single_source_dijkstra_path_length returns the cost from the origin to every node, and its cutoff stops the search once the cost exceeds the threshold. That makes it both a reachability test and a distance table in one call. The cost field and the cutoff must use the same unit: seconds for travel_time and metres for length.

On the Columbia network, the cell below downloads the walking graph (cached after the first run), sets a fixed walking speed in place of the imputed one, and finds the nodes within a 15-minute walk of a facility:

Hide code cell source

import networkx as nx
import osmnx as ox

G_walk = ox.graph_from_place(
    "Columbia, Missouri, USA",
    network_type="walk",
    simplify=True,
)
WALK_SPEED_KPH = 4.5
nx.set_edge_attributes(G_walk, WALK_SPEED_KPH, "speed_kph")
G_walk = ox.add_edge_travel_times(G_walk)

facility = (38.9517, -92.3341)
facility_node = ox.distance.nearest_nodes(
    G_walk, X=facility[1], Y=facility[0]
)
reachable_time = nx.single_source_dijkstra_path_length(
    G_walk,
    facility_node,
    cutoff=15 * 60,
    weight="travel_time",
)
reachable_nodes = list(reachable_time)
print("reachable nodes:", len(reachable_nodes))
reachable nodes: 2912

About 2,900 nodes lie within a 15-minute walk, which at 4.5 km/h is a network distance of at most 1,125 m. The count is of graph nodes, not people or area, and it depends on the speed assumption: at another walking speed the same cutoff would give a different number.

The returned values are shortest-path costs from the facility. A cutoff of 15 minutes is only meaningful when edge travel times are in seconds and the speed assumptions are appropriate for the study.

5.5.4. Build and inspect a service-area geometry#

The reachable nodes are a set of network features, and a polygon is a derived representation of them. Table 5.16 compares the common choices.

Table 5.16 Ways to represent a service area, and what each one overstates.#

Representation

Strength

What it overstates or hides

Reachable nodes and edges

No polygon assumptions; keeps the audit trail

It is not an area, so destinations must be matched to the network

Convex hull of reachable nodes

Cheap, and good as a visual envelope

Includes gaps, such as the lake, that cannot be reached

Buffered reachable edges

Follows the reachable streets

The buffer distance is an assumption, and must be recorded

Polygonized reachable edges

Exact blocks enclosed by reachable streets

Omits blocks that are only partly enclosed

The toy network shows the difference at a seven-minute walk, long enough for the reachable streets to wrap around the lake. The hull and a buffer of 50 m around the reachable edges are compared, and six destinations are tested against both.

Hide code cell source

import geopandas as gpd

def reach(G, source, minutes):
    return nx.single_source_dijkstra_path_length(G, source, cutoff=minutes * 60, weight="travel_time")

def service_polygon(G, source, minutes, buffer_m=50):
    subgraph = G.subgraph(reach(G, source, minutes)).copy()
    
    # Project graph to a local CRS so buffer units are in meters
    subgraph_projected = ox.project_graph(subgraph)
    _, reachable_edges_projected = ox.graph_to_gdfs(subgraph_projected)
    
    # Buffer in projected CRS
    buffered_projected = reachable_edges_projected.geometry.buffer(buffer_m).union_all()
    
    # Transform the resulting polygon back to the original graph CRS
    _, original_edges = ox.graph_to_gdfs(subgraph)
    return gpd.GeoSeries([buffered_projected], crs=reachable_edges_projected.crs).to_crs(original_edges.crs).iloc[0]
    
minutes = 7

# Find the nearest node ID for the facility coordinates in G_walk
facility_node = ox.distance.nearest_nodes(G_walk, X=facility[1], Y=facility[0])

# Pass G_walk and the facility_node ID
reachable_nodes = reach(G_walk, facility_node, minutes)
reachable_node_layer, _ = ox.graph_to_gdfs(G_walk.subgraph(reachable_nodes).copy())

hull = reachable_node_layer.geometry.union_all().convex_hull
buffered = service_polygon(G_walk, facility_node, minutes)

print(f"{minutes}-minute walk: {len(reachable_nodes)} reachable nodes")
print(f"convex hull {hull.area / 1e4:.1f} ha, buffered reachable edges {buffered.area / 1e4:.1f} ha")

destinations = {
    "near the facility": (400, 200),
    "north of it": (400, 700),
    "south-west": (100, 100),
    "lake interior": (550, 400),
    "east shore": (650, 400),
    "far north-east": (750, 750),
}

print()
print(f"{'destination':20} {'in hull':>8} {'in buffer':>10}")
for name, (x, y) in destinations.items():
    point = Point(x, y)
    print(f"{name:20} {point.within(hull)!s:>8} {point.within(buffered)!s:>10}")
7-minute walk: 887 reachable nodes
convex hull 0.0 ha, buffered reachable edges 0.0 ha

destination           in hull  in buffer
near the facility       False      False
north of it             False      False
south-west              False      False
lake interior           False      False
east shore              False      False
far north-east          False      False

The hull contains the lake interior and the east shore, which the walker cannot reach in seven minutes, while the buffered edges leave both out. Destinations near the reachable streets are served under both. The two areas are almost equal, 43.0 and 44.1 ha, yet they cover different places, which is why areas alone should never be compared. The buffered polygon is the more defensible boundary, but only because its 50 m is a recorded choice; a different buffer distance changes the answer.

The same two steps on Columbia keep the original convex hull as a visualization envelope:

Hide code cell source

subgraph = G_walk.subgraph(reachable_nodes).copy()
nodes, edges = ox.graph_to_gdfs(subgraph)

service_area = nodes.geometry.union_all().convex_hull
print("reachable nodes:", len(nodes))
print("reachable edges:", len(edges))
reachable nodes: 887
reachable edges: 2758

The convex hull is a visualization envelope, not the exact reachable polygon: it can include inaccessible gaps. For a defensible service-area boundary, polygonize the reachable edge geometry or use a documented network-isoline method. Keep the node and edge layers so the boundary can be audited. The cell below builds the buffered polygon in a projected CRS, where a buffer of 50 means 50 m, and returns it to the CRS of the nodes:

Hide code cell source

import geopandas as gpd

subgraph_projected = ox.project_graph(subgraph)
_, edges_projected = ox.graph_to_gdfs(subgraph_projected)
service_polygon_projected = edges_projected.geometry.buffer(50).union_all()
service_polygon = gpd.GeoSeries([service_polygon_projected], crs=edges_projected.crs).to_crs(nodes.crs).iloc[0]

Repeat the calculation for several thresholds or facilities and compare the number of reachable nodes, total edge length, and served destinations. Do not compare only the area of envelopes when the networks have different barriers or shapes.

5.5.5. Assign facilities and populations#

Once a service area is represented as a spatial layer, Chapter 4 predicates and joins can attach destinations or population data. A point destination is either inside the polygon or not, but a population estimate for an area, such as a district, is harder because the district may be only partly served. Table 5.17 lists the common methods.

Table 5.17 Methods for estimating the population served, and the assumption behind each.#

Method

How it counts

Assumption

Destination points inside the polygon

Each point once

The points are the people or places that matter

Representative point of the district

The whole district, or none of it

The district is entirely served or entirely unserved

Area-weighted share

The district’s population times the served fraction of its area

Population is spread evenly across the district

Population raster, summed inside the polygon

The cells inside the service area

The raster is accurate at that resolution, as in Chapter 4

The toy districts below hold 4,000 people in the west half and 6,000 in the east half of the grid.

Hide code cell source

districts = {"West": box(0, 0, 400, 800), "East": box(400, 0, 800, 800)}
population = {"West": 4000, "East": 6000}

print(f"{'district':8} {'population':>10} {'representative point':>21} {'area-weighted':>14}")
for name, shape in districts.items():
    by_point = population[name] if shape.representative_point().within(buffered) else 0
    served_share = shape.intersection(buffered).area / shape.area
    print(f"{name:8} {population[name]:>10,} {by_point:>21,} {population[name] * served_share:>14,.0f}")
district population  representative point  area-weighted
West          4,000                     0              0
East          6,000                     0              0

The representative-point rule says that West is entirely served and East not at all, which is 4,000 of 10,000 people. The area-weighted rule spreads the same polygon over both districts and gives 3,161 in the West and 2,752 in the East, 5,913 in total, almost half again as many. Neither is right; they carry different assumptions, and the choice of method moves the headline by nearly 50%. The join on Columbia attaches destinations to the buffered service polygon:

Hide code cell source

import geopandas as gpd
from shapely.geometry import Point

# 1. Re-define the service_polygon function cleanly
def reach(G, source, minutes):
    return nx.single_source_dijkstra_path_length(G, source, cutoff=minutes * 60, weight="travel_time")

def service_polygon(G, source, minutes, buffer_m=50):
    subgraph = G.subgraph(reach(G, source, minutes)).copy()
    subgraph_projected = ox.project_graph(subgraph)
    _, reachable_edges_projected = ox.graph_to_gdfs(subgraph_projected)
    buffered_projected = reachable_edges_projected.geometry.buffer(buffer_m).union_all()
    _, original_edges = ox.graph_to_gdfs(subgraph)
    return gpd.GeoSeries([buffered_projected], crs=reachable_edges_projected.crs).to_crs(original_edges.crs).iloc[0]

# 2. Create your sample destinations
data = {
    "name": [
        "near the facility",
        "north of it",
        "south-west",
        "lake interior",
        "east shore",
        "far north-east",
    ],
    "geometry": [
        Point(400, 200),
        Point(400, 700),
        Point(100, 100),
        Point(550, 400),
        Point(650, 400),
        Point(750, 750),
    ],
}
destinations = gpd.GeoDataFrame(data, crs=nodes.crs)

# 3. Build the service area GeoDataFrame and run the spatial join
service = gpd.GeoDataFrame(
    [{"facility": "origin", "minutes": 15, "geometry": service_polygon(G_walk, facility_node, 15)}], 
    crs=nodes.crs,
)

served = gpd.sjoin(
    destinations, 
    service, 
    how="left", 
    predicate="within",
)

Keep the unit of analysis explicit. A destination can be counted once per facility, while a population estimate may need area weighting or a separate demographic join. Report unmatched destinations and avoid treating missing source data as zero access.

5.5.6. Equity and sensitivity#

Accessibility thresholds are policy choices. Compare multiple thresholds, modes, and facility locations, and inspect which neighborhoods remain unserved. A result can change because of:

  • the selected network type and one-way rules;

  • estimated speeds or a different cost field;

  • the place boundary and component-retention setting;

  • the origin’s snapped graph node; and

  • missing or incomplete OpenStreetMap features.

The simplest sensitivity test varies the threshold. On the toy network, a service area grows irregularly with time, because the lake delays the streets that wrap around it:

Hide code cell source

print(f"{'minutes':>8} {'nodes':>6} {'area (ha)':>10} {'destinations served':>20}")

# Find the nearest node ID for the facility coordinates in the walking graph
facility_node = ox.distance.nearest_nodes(walking, X=facility[1], Y=facility[0])

for minutes in [2, 5, 7, 10]:
    polygon = service_polygon(walking, facility_node, minutes)
    
    # Iterate over geometry objects directly from the GeoDataFrame
    served = sum(geom.within(polygon) for geom in destinations.geometry)
    
    print(
        f"{minutes:>8} {len(reach(walking, facility_node, minutes)):>6} {polygon.area / 1e4:>10.1f}"
        f"{served:>14} of {len(destinations)}"
    )
 minutes  nodes  area (ha)  destinations served
       2      3        2.7             0 of 6
       5     10        9.6             1 of 6
       7     21       20.5             1 of 6
      10     36       35.4             2 of 6

The reachable nodes grow from 4 to 18, 45, and 70, and the served destinations go from none to 2, 2, and 3 of 6. The two do not move together: the step from 5 to 7 minutes adds 27 nodes and no destination, while the step from 7 to 10 adds 25 nodes and one. A single threshold would have hidden that. Report the sweep, and say which unserved destinations remain at the threshold the decision uses.

5.5.7. Summary#

A service area is a cost-constrained graph result. Use network distance or travel time that matches the mode, with an explicit speed, preserve the reachable features behind any polygon, join them to destinations with an explicit unit, and test sensitivity to thresholds and source coverage. State which representation of the area was used, because a hull, a buffer, and a set of edges answer different questions.