Build an Origin-Destination Matrix in PyQGIS

An origin-destination matrix — every origin's network cost to every destination — is the input to catchment allocation, facility siting, accessibility scoring and most of transport planning. It is also the workload that separates people who understand the routing API from people who call the Processing algorithm in a nested loop and wait four hours for a result that should have taken forty seconds.

This recipe belongs to Network Analysis & Routing in PyQGIS. It covers the one-build-many-queries structure, writing the result out in both wide and long form, handling unreachable pairs, and keeping a large matrix affordable.

Where the time goes in a matrixBuilding the graph happens once. A Dijkstra run happens once per origin and fills the cost for every destination at the same time. Reading an individual pair is a single array lookup. So a matrix of fifty origins by two thousand destinations costs one build and fifty runs, not a hundred thousand routings.50 origins × 2,000 destinations = 50 runs, not 100,000build graphonce — seconds to minutestie every point heredijkstra per origin50 runs — milliseconds eachfills every destination at oncearray lookup per pair100,000 — freecosts[destination_vertex]the shape that goes wrongfor origin in origins: for destination in destinations: processing.run(...)100,000 graph builds — hours, for the same answerthe same reasoning applies to service areas, nearest-facility and allocationif the network is fixed, build it once

Prerequisites

  • QGIS 3.34 LTR or newer, everything in one projected CRS.
  • Origin and destination point layers with stable identifiers.
  • The tie mechanics from snapping points to a network.

Tie everything in one build

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

roads = QgsProject.instance().mapLayersByName("roads")[0]
depots = QgsProject.instance().mapLayersByName("depots")[0]
sites = QgsProject.instance().mapLayersByName("sites")[0]

origins = [(f["depot_id"], f.geometry().asPoint()) for f in depots.getFeatures()]
destinations = [(f["site_id"], f.geometry().asPoint()) for f in sites.getFeatures()]

director = QgsVectorLayerDirector(
    roads, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
)
director.addStrategy(
    QgsNetworkSpeedStrategy(roads.fields().indexOf("speed_kph"), 48.0, 1000.0 / 3600.0)
)

points = [p for _, p in origins] + [p for _, p in destinations]
builder = QgsGraphBuilder(roads.sourceCrs(), True, 0.5)
tied = director.makeGraph(builder, points)
graph = builder.graph()

origin_vertices = {
    key: graph.findVertex(tied[i]) for i, (key, _) in enumerate(origins)
}
destination_vertices = {
    key: graph.findVertex(tied[len(origins) + i])
    for i, (key, _) in enumerate(destinations)
}

Breakdown: Concatenating both point lists into one makeGraph call is the whole trick — one build ties everything, and the offset arithmetic (len(origins) + i) recovers which tied point belongs to which layer. Keeping the identifiers alongside is what makes the result joinable back to the source data afterwards; feature ids would not survive a re-export. If either layer is large, this is also where the memory goes, since each tied point splits an edge.

One run per origin

matrix = {}
for origin_key, origin_vertex in origin_vertices.items():
    tree, costs = QgsGraphAnalyzer.dijkstra(graph, origin_vertex, 0)
    matrix[origin_key] = {
        dest_key: (costs[v] if tree[v] != -1 or v == origin_vertex else None)
        for dest_key, v in destination_vertices.items()
    }
    reached = sum(1 for value in matrix[origin_key].values() if value is not None)
    print(f"{origin_key}: {reached}/{len(destination_vertices)} reachable")

Breakdown: The dictionary comprehension reads every destination out of one cost array, so the inner loop is array indexing rather than routing. The reachability test uses tree[v] != -1 rather than a cost threshold, with the origin vertex special-cased because its own tree entry is -1 while its cost is legitimately zero. Storing None for unreachable pairs rather than a sentinel keeps the distinction between "far" and "impossible" intact all the way to the output, which matters enormously when someone later takes a mean.

Printing the reachable count per origin is the diagnostic that catches a fragmented network early — an origin that reaches 12 of 2,000 destinations is on an island, not in a badly connected suburb.

Writing it out

Wide reads well; long joins wellThe wide format puts one row per origin and one column per destination, which is readable in a spreadsheet and awkward past a few dozen columns. The long format puts one row per origin destination pair with the cost as a single column, which is what a database join and most analysis tools expect.Write long unless a person is going to read itwideorigin, site_1, site_2, site_3depot_a, 412, 880, 1304depot_b, 906, 344, 771depot_c, 1120, 690, 288readable; 2,000 columns is notlongorigin, destination, secondsdepot_a, site_1, 412depot_a, site_2, 880depot_a, site_3, 1304joins, filters and aggregates cleanlya blank cell and a zero mean different things — write empty, never 0

import csv

with open("/data/output/od_matrix.csv", "w", newline="") as handle:
    writer = csv.writer(handle)
    writer.writerow(["origin", "destination", "seconds"])
    for origin_key, row in matrix.items():
        for dest_key, seconds in row.items():
            writer.writerow([
                origin_key, dest_key,
                "" if seconds is None else round(seconds, 1),
            ])

Breakdown: Writing an empty string for unreachable pairs rather than 0, -1 or 9999 is the single most important decision in this file. A zero reads as "adjacent" to every tool that consumes it and will silently become the nearest facility for every origin. newline="" is required for correct line endings on Windows and harmless elsewhere. Rounding to one decimal keeps the file small without losing anything meaningful given the accuracy of the underlying speed model.

For the nearest-facility question, which is what most matrices are actually built to answer, collapse as you go:

nearest = {}
for origin_key, row in matrix.items():
    reachable = {k: v for k, v in row.items() if v is not None}
    if reachable:
        best = min(reachable, key=reachable.get)
        nearest[origin_key] = (best, reachable[best])
    else:
        nearest[origin_key] = (None, None)

Breakdown: Filtering out the None values before min is required — comparing None against a float raises in Python 3, which is a good thing here because it forces the unreachable case to be handled rather than silently winning. Keeping both the identifier and the cost lets you write a result layer where each origin is styled by the facility that serves it, which is the map this analysis usually ends up as.

Joining the result back to the map

A CSV is the durable artefact; a layer is what people look at. Adding the nearest-facility result as fields on the origin layer is a short loop.

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

depots.startEditing()
for name, kind in (("nearest_site", QVariant.String), ("seconds", QVariant.Double)):
    if depots.fields().indexOf(name) == -1:
        depots.addAttribute(QgsField(name, kind))
depots.updateFields()

site_index = depots.fields().indexOf("nearest_site")
cost_index = depots.fields().indexOf("seconds")

for feature in depots.getFeatures():
    best, seconds = nearest.get(feature["depot_id"], (None, None))
    depots.changeAttributeValue(feature.id(), site_index, best)
    depots.changeAttributeValue(feature.id(), cost_index, seconds)

depots.commitChanges()

Breakdown: Checking for the field before adding it makes the script safe to re-run, which it will be — matrices get rebuilt with different speed assumptions far more often than anyone expects. updateFields() must come after the additions and before the index lookups, or both indices come back as -1 and every write silently targets the wrong column. Assigning None for an unreachable origin leaves the attribute null rather than zero, which carries the same distinction into the layer that the CSV preserved, and lets a categorized renderer show unserved origins in their own colour rather than lumping them with the fastest.

For the full matrix rather than the collapsed version, a CSV joined to the layer with QgsVectorLayerJoinInfo, or simply loaded as a non-spatial table, is usually better than adding two thousand columns to a point layer.

Keeping a large matrix affordable

The matrix itself is cheap; the graph and the tied points are not. Three levers matter.

Clip the network. A matrix inside one city does not need a national road layer. Clip to the origins and destinations' bounding box plus a buffer wider than the longest journey you expect.

Tie fewer points. Where destinations cluster — a hundred addresses on one street — snapping them all splits that edge a hundred times. Aggregating destinations to a sensible unit before routing, and distributing the result afterwards, is both faster and usually more honest about the accuracy on offer.

Store long, not wide. A wide matrix of 5,000 × 5,000 floats is 200 MB in memory as Python objects and far more as a CSV of mostly-repeated identifiers. Streaming rows to the file as each origin finishes, instead of accumulating a dictionary, keeps memory flat regardless of size.

QGIS version compatibility

Everything here uses qgis.analysis classes that have been stable since QGIS 3.0, with the fromVertex()/toVertex() rename in 3.24 the only change of note — and this recipe does not walk edges, so it is unaffected. QgsGraphAnalyzer.dijkstra has had the same signature throughout.

Troubleshooting

  • The script takes hours. A Processing call inside the loop, rebuilding the graph per pair.
  • Every cost is None. The destinations were tied against a different build than the one being queried.
  • Costs are asymmetric on a two-way network. They should not be — check for a direction field whose default is forward rather than both.
  • A mean travel time is implausibly low. Unreachable pairs were written as 0 and are dragging it down.
  • Memory climbs until the process dies. The whole matrix is being accumulated in a dictionary; stream to the file instead.
  • Identifiers do not match the source layer afterwards. Feature ids were used instead of a stable business key.

Conclusion

Tie every origin and destination in one graph build, run Dijkstra once per origin, read every destination out of the resulting cost array, and write the result long with unreachable pairs left empty. That structure turns a job people expect to take hours into one that takes under a minute on a city network, and the discipline about None is what keeps the answer trustworthy afterwards.

Frequently Asked Questions

Should I run Dijkstra from origins or destinations? From whichever side has fewer points, since the cost is one run per seed. On a symmetric two-way network the matrix is the transpose either way; on a one-way network, running from destinations gives you the inward matrix, which is a different question.

Can I parallelise the origin loop? Not safely across threads with a shared graph in the Python bindings. Splitting origins across separate processes, each building its own graph, works and is worth it only when the origins are many and the network is small.

How do I include the snap offsets in the costs? Add each origin's and destination's offset, converted to time at a walking speed, to the cost. Do it explicitly in the output rather than silently, so the two components stay separable.

What if I need distance and time together? Add both strategies to the director and run Dijkstra twice per origin with criterion 0 and 1. That is two runs per origin against one build, still nothing like the cost of rebuilding.