Find the Nearest Facility by Network Distance in PyQGIS

"Which fire station is closest to this address?" sounds like a nearest-neighbour question, and in straight-line terms it is. On a road network the answer is often different: a river with one bridge, a motorway without a junction, a one-way system — any of them can make the geometrically nearest station the slower one to arrive from. Network nearest-facility analysis assigns each location to the facility with the lowest travel cost along the roads, which is what planning, emergency response and service allocation actually need.

This recipe belongs to Network Analysis & Routing. It builds a graph once, runs one shortest-path tree per facility rather than one route per location, keeps the best facility for every location, compares results with straight-line assignment, and handles unreachable locations and facility capacity.

Straight-line nearest versus network nearestAn address beside a river. Station A is 600 metres away in a straight line, but across the river with the nearest bridge two kilometres upstream, so its road distance is 4.1 kilometres. Station B is 1.4 kilometres away in a straight line and 1.8 kilometres by road on the same bank. Network analysis assigns the address to B.Nearest by crow, or by road?riverbridgeaddressA: 0.6 km straight, 4.1 km roadB: 1.4 km straight, 1.8 km road

Prerequisites

Build the graph with every point tied in

Facilities and locations must become graph vertices. Passing all of them to the director as tie points snaps each onto the network once.

from qgis.core import QgsProject, QgsVectorLayer
from qgis.analysis import (QgsVectorLayerDirector, QgsNetworkDistanceStrategy,
                           QgsGraphBuilder, QgsGraphAnalyzer)

roads = QgsProject.instance().mapLayersByName("roads")[0]
stations = QgsProject.instance().mapLayersByName("fire_stations")[0]
addresses = QgsProject.instance().mapLayersByName("addresses")[0]

director = QgsVectorLayerDirector(roads, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth)
director.addStrategy(QgsNetworkDistanceStrategy())
builder = QgsGraphBuilder(roads.crs(), True, 1.0)

st_pts = [(f["station_id"], f.geometry().asPoint()) for f in stations.getFeatures()]
ad_pts = [(f.id(), f.geometry().asPoint()) for f in addresses.getFeatures()]
tied = director.makeGraph(builder, [p for _, p in st_pts] + [p for _, p in ad_pts])
graph = builder.graph()

st_vertices = [(sid, graph.findVertex(tied[i])) for i, (sid, _) in enumerate(st_pts)]
offset = len(st_pts)
ad_vertices = [(fid, graph.findVertex(tied[offset + i])) for i, (fid, _) in enumerate(ad_pts)]
print(graph.vertexCount(), "vertices;", len(st_vertices), "stations;", len(ad_vertices), "addresses")

Breakdown: makeGraph returns the tied (snapped) positions in the same order as the input points, so the first entries belong to stations and the rest to addresses; keeping that order straight is the main bookkeeping in this workflow. The topology tolerance of one metre merges nearly coincident road ends. Distance strategy makes edge cost equal to length; for travel time, use a speed strategy as in adding speed and travel cost to a network. Tying every address into the graph is fine for thousands of points; for hundreds of thousands, snap addresses to the nearest graph vertex instead.

Run one tree per facility

Routing from each address to each station would mean addresses × stations shortest paths. Dijkstra from a station computes the cost to every vertex in one run, so one run per station covers all addresses.

One shortest-path tree per facilityFor each of a few stations, Dijkstra computes the cost from that station to every vertex of the graph in one run. The address vertices' costs are read from each tree, and for every address the station with the lowest cost is kept. Ten stations and fifty thousand addresses need ten runs, not half a million routes.Facilities × 1 run, not facilities × locationsstation sDijkstracost to all verticesread address costscost[v] for eachkeep minimumbest stationbest cost

import math

best = {fid: (None, math.inf) for fid, _ in ad_vertices}
for sid, sv in st_vertices:
    tree, cost = QgsGraphAnalyzer.dijkstra(graph, sv, 0)
    for fid, av in ad_vertices:
        c = cost[av]
        if tree[av] != -1 or av == sv:          # reachable
            if c < best[fid][1]:
                best[fid] = (sid, c)

unreached = [fid for fid, (sid, c) in best.items() if sid is None]
print(len(unreached), "addresses not reachable from any station")

Breakdown: dijkstra returns two lists indexed by vertex: the incoming edge of the shortest-path tree and the cost from the start. A vertex is reachable if it has an incoming edge (or is the start itself); unreachable vertices keep an infinite or meaningless cost and must be excluded. Updating the best station per address as each tree is computed keeps memory low — only one tree exists at a time. Ten stations and fifty thousand addresses take ten Dijkstra runs, typically seconds.

Write the assignment to a layer

The result belongs on the addresses: the assigned station and the network distance, ready to map, count and report.

from qgis.core import QgsField, edit
from qgis.PyQt.QtCore import QVariant

prov = addresses.dataProvider()
for name, qtype in (("nearest_station", QVariant.String), ("road_km", QVariant.Double)):
    if addresses.fields().indexOf(name) < 0:
        prov.addAttributes([QgsField(name, qtype)])
addresses.updateFields()
i_st, i_km = addresses.fields().indexOf("nearest_station"), addresses.fields().indexOf("road_km")

changes = {fid: {i_st: sid, i_km: (round(c / 1000, 3) if sid else None)}
           for fid, (sid, c) in best.items()}
prov.changeAttributeValues(changes)
print("assigned", sum(1 for s, _ in best.values() if s), "addresses")

Breakdown: One provider call writes both fields for every address. Unreachable addresses get NULL rather than a misleading station. Styling addresses by nearest_station with a categorized renderer shows each station's service area as a coloured point pattern; dissolving Voronoi-like regions from those points, or calculating service areas per station, turns it into polygons.

Compare with straight-line assignment

The interesting addresses are those where the network answer differs from the straight-line answer — they show where barriers shape access.

from qgis.core import QgsSpatialIndex

st_index = QgsSpatialIndex()
st_geoms = {}
for f in stations.getFeatures():
    st_index.addFeature(f)
    st_geoms[f.id()] = f["station_id"]

differs = 0
for f in addresses.getFeatures():
    nearest_fid = st_index.nearestNeighbor(f.geometry().asPoint(), 1)[0]
    if st_geoms[nearest_fid] != best[f.id()][0]:
        differs += 1
print(f"{differs} addresses ({differs / addresses.featureCount():.1%}) are served by a "
      "different station than the straight-line nearest")

Breakdown: A spatial index finds the straight-line nearest station quickly; comparing it with the network result counts addresses where geometry misleads. Mapping them usually reveals rivers, railways and motorways as clear boundaries. The share is a useful headline for anyone who suggests that straight-line distance is "good enough".

Report coverage against a standard

Service planning is usually judged against a standard: 90 % of addresses within 8 minutes of a fire station, every pupil within 3 km of a primary school. With the best cost per address computed, coverage is a count.

STANDARD_KM = 3.0
within = [fid for fid, (sid, c) in best.items() if sid and c / 1000 <= STANDARD_KM]
share = len(within) / max(len(best), 1)
print(f"{share:.1%} of addresses within {STANDARD_KM} km road distance of a station")

by_station = {}
for fid, (sid, c) in best.items():
    if sid:
        by_station.setdefault(sid, []).append(c / 1000)
for sid, d in sorted(by_station.items()):
    d.sort()
    print(f"{sid:<10} {len(d):>6} addresses, median {d[len(d) // 2]:.1f} km, "
          f"90th pct {d[int(len(d) * 0.9)]:.1f} km")

Breakdown: The share of addresses within the standard is the headline figure; per-station counts and distance percentiles show which stations carry most of the load and which serve a long tail of remote addresses. Mapping the addresses outside the standard identifies where a new facility would help most — they are the natural candidates to feed into a location-allocation exercise. For time-based standards, build the graph with a speed strategy so costs are already in seconds.

Respect facility capacity

Nearest assignment ignores capacity: a school with 300 places may be nearest for 600 pupils. A simple greedy pass assigns locations in order of how much they lose by not getting their first choice.

Capacity-aware assignmentEach location has a cost to every facility. Locations are sorted by regret, the difference between their best and second-best costs. In that order, each is assigned to its cheapest facility with remaining capacity. Locations that would lose most by moving are placed first, giving a reasonable assignment without an optimisation solver.Assign the most constrained firstcost matrixlocation × facilitysort by regret2nd best − bestassigncheapest with room

costs = {fid: {} for fid, _ in ad_vertices}
for sid, sv in st_vertices:
    tree, cost = QgsGraphAnalyzer.dijkstra(graph, sv, 0)
    for fid, av in ad_vertices:
        if tree[av] != -1 or av == sv:
            costs[fid][sid] = cost[av]

capacity = {f["station_id"]: f["capacity"] for f in stations.getFeatures()}
def regret(fid):
    c = sorted(costs[fid].values())
    return (c[1] - c[0]) if len(c) > 1 else math.inf
assigned = {}
for fid in sorted(costs, key=regret, reverse=True):
    for sid, c in sorted(costs[fid].items(), key=lambda kv: kv[1]):
        if capacity[sid] > 0:
            assigned[fid] = sid
            capacity[sid] -= 1
            break
print(len(assigned), "assigned within capacity")

Breakdown: Keeping the full cost row per location — one cost per facility — allows second choices. Ordering by regret places first the locations that would suffer most from not getting their nearest facility, a classic heuristic that gives sensible results quickly. It is not guaranteed optimal; for formal planning, export the cost matrix to an optimisation library. The origin–destination matrix recipe builds the same matrix in table form.

QGIS version compatibility

QgsVectorLayerDirector, QgsGraphBuilder and QgsGraphAnalyzer.dijkstra are available in qgis.analysis on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. On QGIS 4, QgsVectorLayerDirector.DirectionBoth is QgsVectorLayerDirector.Direction.DirectionBoth.

Troubleshooting

  • Many addresses are unreachable. The network has disconnected pieces; check connectivity and snap undershoots.
  • Every address goes to one station. Station vertices were mixed up with address vertices; check the order of tied points.
  • Runs are slow. Too many tie points; snap addresses to existing vertices instead of tying each in.
  • Distances look too short. Cost is in map units of a geographic CRS; reproject to metres.

Conclusion

Tie facilities and locations into one graph, run one Dijkstra tree per facility and keep each location's minimum cost, write the station and distance back in one call, compare with straight-line assignment to reveal barriers, and use a regret-ordered greedy pass when facilities have limited capacity.

Frequently Asked Questions

Can I use travel time instead of distance? Yes, with a speed strategy on the director; costs then come out in seconds or hours.

What about one-way streets? Set the director's direction field and values so the graph respects one-way restrictions.

Is there a Processing algorithm for this?native:shortestpathpointtolayer and service-area algorithms cover related cases; a script is clearer for many-to-many assignment.

How many facilities can this handle? Each facility is one Dijkstra run; hundreds are fine.