How to Build a Distance Matrix Between Two Layers in Python

On this page

Problem statement

"How far is every home from every bus stop?" is a natural question with an expensive answer. For Leeds' 2,607 neighbourhoods against West Yorkshire's 16,450 bus stops it is 42.9 million distances, 343 MB as a float64 array (see The data used on this page). Do the same for all 188,880 output areas in England and Wales against that many stops and it is 3.1 billion distances — 25 GB.

The loop most people write first is worse than the memory:

distances = []
for home in homes.geometry:
    distances.append(stops.distance(home).values)
matrix = np.vstack(distances)

It builds the same matrix, one row at a time, through the Python interpreter — and it almost never needed the full matrix in the first place. Nearly every real question is one of three:

the nearest one — which stop is closest, and how far
the k nearest — the three closest, for a choice model
everything within a radius — all stops within 800 m

Each of those has an exact answer that is orders of magnitude cheaper than the full matrix, because a spatial index lets you skip pairs that cannot possibly qualify.

Quick answer

Pick the tool from the question, not from the data:

import numpy as np
from scipy.spatial import cKDTree

a = homes.to_crs("EPSG:27700")            # both layers, same projected CRS
b = stops.to_crs("EPSG:27700")

xa = np.column_stack([a.geometry.x, a.geometry.y])
xb = np.column_stack([b.geometry.x, b.geometry.y])

tree = cKDTree(xb)
distance, index = tree.query(xa, k=1)       # nearest stop to each home
a["nearest_school"] = b.iloc[index].index.values
a["distance_m"] = distance

Question	Tool	Cost for 2,607 × 16,450
nearest	`cKDTree.query(k=1)` or `gpd.sjoin_nearest`	0.003 s / 0.007 s
k nearest	`cKDTree.query(k=3)`	0.003 s
within a radius	`cKDTree.sparse_distance_matrix`	0.004 s, 37,623 pairs within 500 m
genuinely all pairs	`scipy.spatial.distance.cdist`	0.14 s and 343 MB here; grows as the product

Only the last row needs the full matrix, and it is the row you almost never actually want.

What you get — from 2. For the nearest neighbour, use sjoin_nearest or a KD-tree

       id  index_right  distance_m
E00056750         8162   40.853396
E00056751         5688  275.261694
E00056752         8159  116.077560

One CRS, then nearest, k nearest, radius or all pairs

1. Put both layers in the same projected CRS

if a.crs != b.crs:
    b = b.to_crs(a.crs)
if a.crs.is_geographic:
    raise ValueError("distances in degrees are not distances")

A KD-tree computes Euclidean distance on whatever numbers you give it. On degrees that is neither metres nor consistent — 0.01 degrees is 1.1 km north-south and 660 m east-west in Leeds. See how to measure distance accurately.

2. For the nearest neighbour, use `sjoin_nearest` or a KD-tree

GeoPandas has this built in, and it keeps the attributes:

joined = gpd.sjoin_nearest(a, b, how="left", distance_col="distance_m")
print(joined[["id", "index_right", "distance_m"]].head(3).to_string(index=False))

       id  index_right  distance_m
E00056750         8162   40.853396
E00056751         5688  275.261694
E00056752         8159  116.077560

The KD-tree version is faster on large inputs and gives you the index directly:

tree = cKDTree(xb)
distance, index = tree.query(xa, k=1)

They agree exactly — on the Leeds homes and stops the maximum difference between the two was 0.00e+00 m (measured under Why a KD-tree is so much faster). Use sjoin_nearest when you want the joined attributes, the tree when you want speed and control.

3. For k nearest, ask for k

distance, index = tree.query(xa, k=3)
print(distance[:2].round(1))
print(index[:2])

[[ 40.9  51.2  81.2]
 [275.3 285.8 377.1]]
[[ 8162  8163 16030]
 [ 5688  5687 11913]]

Both arrays come back with shape (len(xa), k), sorted nearest first. This is the shape a choice model or an accessibility index wants, and it costs barely more than k=1.

4. For a radius, build a sparse matrix

tree_a, tree_b = cKDTree(xa), cKDTree(xb)
sparse = tree_a.sparse_distance_matrix(tree_b, max_distance=500, output_type="coo_matrix")
print(f"{sparse.nnz:,} pairs within 500 m, of {len(xa) * len(xb):,} possible")

37,623 pairs within 500 m, of 42,885,150 possible

Less than a tenth of one percent of the pairs. The sparse matrix holds only those, and it took 0.004 seconds — against 343 MB for the dense equivalent here, and a MemoryError at national scale.

5. Only build the dense matrix when you truly need every pair

from scipy.spatial.distance import cdist

if len(xa) * len(xb) > 50_000_000:
    raise MemoryError(f"{len(xa) * len(xb):,} pairs = {len(xa) * len(xb) * 8 / 1e9:.1f} GB")
matrix = cdist(xa, xb)

cdist is fully vectorised and fast — 42.9 million pairs in 0.14 seconds here. The limit is memory, not time, so guard on the pair count rather than waiting for the crash.

The matrix grows as the product of the two sizes. The KD-tree answered the middle case in 0.003 s without building any of it.

Leeds output areas to their nearest transport stops

The data used on this page

Every example runs on real data: the 2,607 Census 2021 output areas of Leeds as homes — each one a neighbourhood of about 120 households, located at the population-weighted centroid ONS publishes — and the Department for Transport's NaPTAN register of public transport stops for West Yorkshire: the active bus stops as stops, and the railway station entrances as stations. Stops outside Leeds are kept on purpose: the nearest stop to a home near the boundary may be over it.

import io

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Output_Areas_2021_EW_BGC_V2/FeatureServer/0/query")
rows, offset = [], 0
while True:
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "OA21CD,BNG_E,BNG_N",
                                     "returnGeometry": "false", "f": "json", "resultOffset": offset}, timeout=120).json()
    rows += [f["attributes"] for f in page["features"]]
    if not page.get("exceededTransferLimit"):
        break
    offset += len(page["features"])
oa = pd.DataFrame(rows).sort_values("OA21CD").reset_index(drop=True)
homes = gpd.GeoDataFrame({"id": oa["OA21CD"]}, geometry=gpd.points_from_xy(oa["BNG_E"], oa["BNG_N"]), crs=27700)

# NaPTAN area 450 as saved on 9 October 2026 (the register changes weekly); live:
# naptan.api.dft.gov.uk/v1/access-nodes?dataFormat=csv&atcoAreaCodes=450
naptan = pd.read_csv(io.StringIO(requests.get("https://data.spatialworkflow.io/naptan-area-450-2026-10-09.csv",
                                              timeout=300).text), low_memory=False)
naptan = naptan[naptan["Status"] == "active"]
as_points = lambda df: gpd.GeoDataFrame(df[["ATCOCode", "CommonName"]].reset_index(drop=True), crs=27700,
                                        geometry=gpd.points_from_xy(df["Easting"], df["Northing"]))
stops = as_points(naptan[naptan["StopType"].isin(["BCT", "BCS"])])      # on-street bus stops and bus station bays
stations = as_points(naptan[naptan["StopType"] == "RSE"])              # railway station entrances
print(f"{len(homes):,} homes (output areas), {len(stops):,} bus stops, {len(stations):,} station entrances")

2,607 homes (output areas), 16,450 bus stops, 155 station entrances

Example 1 — one function that picks the right method

import geopandas as gpd
import numpy as np
import pandas as pd
from scipy.spatial import cKDTree


def _coords(gdf):
    if not (gdf.geometry.geom_type == "Point").all():
        gdf = gdf.copy()
        gdf["geometry"] = gdf.geometry.representative_point()
    return np.column_stack([gdf.geometry.x, gdf.geometry.y])


def distances(source, target, *, k=None, radius=None, max_dense_pairs=50_000_000):
    """Nearest (default), k-nearest, or all pairs within a radius."""
    if source.crs != target.crs:
        target = target.to_crs(source.crs)
    if source.crs is None or source.crs.is_geographic:
        raise ValueError(f"{source.crs} is geographic — distances would be in degrees")

    xa, xb = _coords(source), _coords(target)
    tree = cKDTree(xb)

    if radius is not None:
        pairs = cKDTree(xa).sparse_distance_matrix(tree, radius, output_type="coo_matrix")
        print(f"{pairs.nnz:,} pairs within {radius} of {len(xa) * len(xb):,} possible "
              f"({pairs.nnz / (len(xa) * len(xb)):.3%})")
        return pd.DataFrame({
            "source": source.index.to_numpy()[pairs.row],
            "target": target.index.to_numpy()[pairs.col],
            "distance": pairs.data,
        }).sort_values(["source", "distance"], ignore_index=True)

    k = k or 1
    distance, index = tree.query(xa, k=k)
    if k == 1:
        distance, index = distance[:, None], index[:, None]

    out = source.copy()
    for j in range(k):
        out[f"target_{j + 1}"] = target.index.to_numpy()[index[:, j]]
        out[f"distance_{j + 1}"] = distance[:, j]
    print(f"{len(out):,} sources, {k} nearest each; "
          f"median nearest {np.median(distance[:, 0]):,.0f}")
    return out


nearest = distances(homes, stops)
three = distances(homes, stops, k=3)
within = distances(homes, stops, radius=500)
print(within.head(3).to_string(index=False))

2,607 sources, 1 nearest each; median nearest 142
2,607 sources, 3 nearest each; median nearest 142
37,623 pairs within 500 of 42,885,150 possible (0.088%)
 source  target  distance
      0    8162 40.853396
      0    8163 51.195703
      0   16030 81.154174

representative_point() rather than centroid for non-point inputs: a centroid can fall outside a concave polygon, which for a distance-to-nearest question puts the reference in the wrong place entirely.

Example 2 — an accessibility measure from the k-nearest result

def accessibility(source, target, *, k=3, decay=800):
    """Sum of distance-decayed access to the k nearest opportunities."""
    result = distances(source, target, k=k)
    columns = [f"distance_{j + 1}" for j in range(k)]
    d = result[columns].to_numpy()

    result["access"] = np.exp(-d / decay).sum(axis=1)
    result["nearest_m"] = d[:, 0]

    print(result["access"].describe()[["min", "50%", "max"]].round(2).to_string())
    return result


scored = accessibility(homes, stops, k=3, decay=800)
worst = scored.nlargest(3, "nearest_m")[["nearest_m", "access"]]
print(worst.round(1).to_string())

2,607 sources, 3 nearest each; median nearest 142
min    0.30
50%    2.41
max    2.88
      nearest_m  access
24       1801.2     0.3
30       1777.8     0.3
2078     1251.2     0.6

The exponential decay is the modelling choice, and decay=800 says access halves roughly every 550 m. Make it a named parameter and state it in the output — it is exactly the kind of number that silently decides a result.

Note that max is 2.88, not 3: even a neighbourhood with three stops close by does not score a full 3, because none of them is at zero distance. At the other end, the least accessible output area is 1.8 km from its nearest stop. If your index needs to be interpretable, normalise it.

Example 3 — a genuine all-pairs matrix, guarded

from scipy.spatial.distance import cdist


def full_matrix(source, target, *, max_gb=2.0):
    if source.crs != target.crs:
        target = target.to_crs(source.crs)

    xa, xb = _coords(source), _coords(target)
    gb = len(xa) * len(xb) * 8 / 1e9
    if gb > max_gb:
        raise MemoryError(
            f"{len(xa):,} x {len(xb):,} = {len(xa) * len(xb):,} pairs ({gb:.1f} GB). "
            f"Use k-nearest or a radius instead, or raise max_gb deliberately."
        )

    matrix = cdist(xa, xb)
    print(f"{matrix.shape} matrix, {matrix.nbytes / 1e6:.0f} MB, "
          f"min {matrix.min():.0f} max {matrix.max():,.0f}")
    return pd.DataFrame(matrix, index=source.index, columns=target.index)


small = full_matrix(stations, homes)
print(small.iloc[:3, :4].round(0).to_string())

(155, 2607) matrix, 3 MB, min 33 max 57,134
         0        1        2        3
0  12288.0  13066.0  12648.0  11949.0
1  12178.0  12956.0  12538.0  11841.0
2  12160.0  12941.0  12522.0  11821.0

155 station entrances by 2,607 homes — 3 MB — is a perfectly reasonable dense matrix, and it is the shape an optimiser wants. The guard exists so that the day someone passes every output area in the country, they get a message naming the alternative rather than a killed process.

Python loops, KD-trees and the radius question

Why the naive loop is slow twice over

for home in homes.geometry:
    stops.distance(home)

Two costs compound. Each .distance() call is a vectorised operation over all of stops — fine — but the loop around it runs in Python, once per home, so you pay interpreter overhead once per home — 2,607 times here. And it computes every pair regardless: no index, no pruning.

Measured below on the Leeds data, the loop took 1.07 s against 0.003 s for the KD-tree — over 300 times slower on a modest problem, and the gap widens with size because the tree prunes and the loop does not.

Why a KD-tree is so much faster

A KD-tree recursively splits the points along alternating axes. To find the nearest neighbour of a query point, it descends to the leaf containing it, then backtracks — but only into branches whose bounding region is closer than the best distance found so far.

For well-distributed 2D points the vast majority of branches are eliminated, giving roughly O(log n) per query rather than O(n). Building the tree is O(n log n) and happens once.

The measured effect on the 2,607 homes and 16,450 stops, best of three:

import time
from scipy.spatial.distance import cdist


def best(fn, repeat=3):
    times = []
    for _ in range(repeat):
        t = time.perf_counter()
        out = fn()
        times.append(time.perf_counter() - t)
    return min(times), out


xa = np.column_stack([homes.geometry.x, homes.geometry.y])
xb = np.column_stack([stops.geometry.x, stops.geometry.y])
t_loop, _ = best(lambda: [stops.distance(h) for h in homes.geometry], repeat=1)
t_full, full = best(lambda: cdist(xa, xb))
t_tree, (d1, i1) = best(lambda: cKDTree(xb).query(xa, k=1))
t_k3, _ = best(lambda: cKDTree(xb).query(xa, k=3))
t_sparse, sp = best(lambda: cKDTree(xa).sparse_distance_matrix(cKDTree(xb), 500, output_type="coo_matrix"))
t_join, joined = best(lambda: gpd.sjoin_nearest(homes, stops, distance_col="d"))
agree = np.abs(joined.groupby(level=0)["d"].first().to_numpy() - d1).max()
print(f" python loop    {t_loop:8.3f} s")
print(f" cdist, all     {t_full:8.3f} s   {full.nbytes / 1e6:,.0f} MB")
print(f" kdtree k=1     {t_tree:8.3f} s")
print(f" kdtree k=3     {t_k3:8.3f} s")
print(f" sparse r<500   {t_sparse:8.3f} s   {sp.nnz:,} pairs")
print(f" sjoin_nearest  {t_join:8.3f} s   max difference from the tree {agree:.2e} m")

 python loop       0.931 s
 cdist, all        0.097 s   343 MB
 kdtree k=1        0.003 s
 kdtree k=3        0.003 s
 sparse r<500      0.004 s   37,623 pairs
 sjoin_nearest     0.006 s   max difference from the tree 0.00e+00 m

Why the radius question is the important one

Most real questions have a threshold in them — a walking distance, a service radius, a catchment. Once a threshold exists, the great majority of pairs are irrelevant, and a sparse structure holds only the relevant ones.

In the example above, 500 m kept 0.088% of the pairs. The dense matrix would have spent 99.9% of its 343 MB storing distances nobody was going to look at.

sparse_distance_matrix returns a COO matrix — parallel arrays of row, column and value — which converts directly into the long-format DataFrame that a groupby or a join wants.

A spatial index computes the second number without ever building the first.

Why straight-line distance may be the wrong answer entirely

Everything here measures Euclidean distance. If a river, a railway or a motorway separates two points, the straight line is not the trip anyone makes.

For questions about travel — accessibility, catchments, service areas — you want network distance along the street graph, which is a different computation with different tools. Straight-line distance is a reasonable proxy in dense, well-connected urban areas and a poor one across any barrier.

Use Euclidean when the question is genuinely about proximity (nearest sensor, spatial weights, clustering) and network distance when it is about travel.

Euclidean only, array shapes and ties

cKDTree uses Euclidean distance only. For lat/lon, use sklearn.neighbors.BallTree with metric="haversine" and radians.
query(k=1) returns 1-D arrays; k>1 returns 2-D. Handle both, or your indexing silently changes shape.
Unmatched queries return inf and n (one past the last index) when distance_upper_bound is set. Check for them before using the index.
sjoin_nearest duplicates rows on ties. Two targets at exactly the same distance both match. Deduplicate if you need one row per source.
Polygons need a representative point. representative_point() is guaranteed inside the shape; centroid is not.
Sparse matrices store explicit zeros for coincident points. A distance of exactly 0 is a real value in COO format and easy to lose when converting.
Memory is the product, not the sum. 188,880 × 16,450 is 25 GB whichever way round you put them.
Guard on the pair count, not the row count. 1,000,000 × 200,000 would be 1,600 GB.

Internal links

Nearest-neighbour joins explained — ties, search radius and the join semantics
How to measure distance accurately in Python — projected, geodesic, and when each is right
Network distance vs straight-line distance explained — when Euclidean is the wrong question
How to use the spatial index in GeoPandas — the same pruning idea inside GeoPandas
How to cluster points by location with DBSCAN in Python — a consumer of the same neighbour queries
Choose a projected CRS for your area — before measuring anything
How to calculate area and distance in GeoPandas correctly — the GeoPandas-level view
How to measure network distance for many origin-destination pairs — the routed equivalent

FAQ

How do I compute distances between every pair of points?

Usually you should not. For every output area in England and Wales against 16,450 stops that is 25 GB. Ask instead for the nearest, the k nearest, or everything within a radius — all of which are far cheaper and exact.

What is the fastest way to find the nearest feature?

scipy.spatial.cKDTree.query(k=1), or gpd.sjoin_nearest if you want the target's attributes joined on. They give identical distances.

How do I get all pairs within a distance?

cKDTree.sparse_distance_matrix(other_tree, max_distance). It returns only the qualifying pairs — 0.088% of them in the example here.

Does this work on latitude and longitude?

Not with cKDTree, which is Euclidean. Project to metres, or use BallTree with metric="haversine" and coordinates in radians.

What about polygons rather than points?

Use representative_point() for a distance-to-nearest question. If you need true polygon-to-polygon distance, gpd.sjoin_nearest computes edge-to-edge distance properly.

Why is my loop over `.distance()` so slow?

It pays Python interpreter overhead once per row and computes every pair with no pruning. A KD-tree does neither.

Should I use straight-line or network distance?

Euclidean for proximity questions — nearest sensor, spatial weights, clustering. Network distance for anything about travel, especially where rivers, railways or motorways form barriers.