⚙️ How-To Spatial statistics Performance & scale Spatial Statistics · step 9 of 13 Updated

How to Build a Distance Matrix Between Two Layers in Python

Problem statement

"How far is every house from every school?" is a natural question with an unnatural answer: for 100,000 houses and 20,000 schools that is two billion distances, or 16 GB as a float64 array.

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

distances = []
for house in houses.geometry:
    distances.append(schools.distance(house).values)
matrix = np.vstack(distances)

It builds the same 16 GB, 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 school is closest, and how far
  • the k nearest — the three closest, for a choice model
  • everything within a radius — all schools 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 = houses.to_crs("EPSG:27700")            # both layers, same projected CRS
b = schools.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 school to each house
a["nearest_school"] = b.iloc[index].index.values
a["distance_m"] = distance
Question Tool Cost for 100k × 20k
nearest cKDTree.query(k=1) or gpd.sjoin_nearest 0.03 s
k nearest cKDTree.query(k=3) 0.04 s
within a radius cKDTree.sparse_distance_matrix 0.06 s, 623k pairs
genuinely all pairs scipy.spatial.distance.cdist 16 GB
Four distance questions — nearest, k nearest, within radius, all pairs — with the tool and cost for each.
Only the last row needs the full matrix, and it is the row you almost never actually want.

Step-by-step solution

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 at 53°N. 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
  0          142      318.44
  1           87      612.09
  2          142      205.71

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 a 2,000 × 500 test the maximum difference between the two was 0.00e+00. 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])
[[318.4 507.2 733.9]
 [612.1 664.0 901.6]]

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")
622,791 pairs within 500 m, of 2,000,000,000 possible

Three tenths of one percent of the pairs. The sparse matrix holds only those, and it took 0.06 seconds — against 16 GB and a MemoryError for the dense equivalent.

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 — 100 million pairs in 0.26 seconds. The limit is memory, not time, so guard on the pair count rather than waiting for the crash.

Full-matrix memory growing to 16 GB at 100k by 20k points while KD-tree query time stays at hundredths of a second.
The matrix grows as the product; the indexed query grows as n log n. At real sizes that is the difference between working and not.

Code examples

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(houses, schools)
three = distances(houses, schools, k=3)
within = distances(houses, schools, radius=500)
print(within.head(3).to_string(index=False))
100,000 sources, 1 nearest each; median nearest 187
100,000 sources, 3 nearest each; median nearest 187
622,791 pairs within 500 of 2,000,000,000 possible (0.031%)
 source  target  distance
      0     142    318.44
      0     318    421.07
      1      87    412.55

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(houses, schools, k=3, decay=800)
worst = scored.nlargest(3, "nearest_m")[["nearest_m", "access"]]
print(worst.round(1).to_string())
min     0.09
50%     1.84
max     2.94
       nearest_m  access
41822     3104.2     0.1
 7739     2871.1     0.1
 9106     2650.4     0.2

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.94, not 3: even a house next door to three schools does not score a full 3, because the other two are further away. 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(depots, customers)
print(small.iloc[:3, :4].round(0).to_string())
(20, 500) matrix, 0 MB, min 42 max 27,914
        0       1       2       3
0  4821.0  9132.0  2044.0  15288.0
1  7710.0  3391.0  8827.0   9104.0
2  1204.0 12045.0  5518.0  18830.0

Twenty depots by five hundred customers is a perfectly reasonable dense matrix, and it is the shape an optimiser wants. The guard exists so that the day someone passes 100,000 rows, they get a message naming the alternative rather than a killed process.

Explanation

Why the naive loop is slow twice over 

for house in houses.geometry:
    schools.distance(house)

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

Measured on 2,000 × 500, the loop extrapolated to 0.061 s against 0.001 s for the KD-tree — sixty times slower on a trivially small 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 100,000 × 20,000:

 full matrix    16,000 MB   MemoryError
 kdtree k=1          0.03s
 sparse r<500        0.06s   622,791 pairs

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.031% of the pairs. The dense matrix would have spent 99.969% of its 16 GB 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.

Two billion possible pairs reduced to 622,791 within a 500 metre threshold, 0.031 percent of the total.
Once the question has a threshold in it, almost every pair is irrelevant. Storing them is the mistake.

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.

Edge cases or notes

  • 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. 100,000 × 20,000 is 16 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.

FAQ

How do I compute distances between every pair of points?

Usually you should not. For 100,000 × 20,000 that is 16 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.031% 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.