Spatial Clustering for Delivery Zone Design

Everybody reaches for k-means first, and on delivery stops it produces zones with straight edges that cut across a river. DBSCAN produces shapes that follow where the work actually is — but only if epsilon comes from the data rather than from guessing.

Why density, not centroids

k-means partitions space by distance to a centre, which means every point belongs to a cluster and the boundaries are the perpendicular bisectors between centres — a Voronoi diagram. For delivery stops that is the wrong model twice over. Real stop distributions are not blobs: they are dense along high streets, dense in apartment complexes, and empty across industrial land and water. And a delivery zone that contains three stops thirty kilometres from the rest of it is not improved by the fact that they were technically closest to that centroid.

DBSCAN, from Ester, Kriegel, Sander and Xu at KDD 1996, asks a different question: which points are in a region of high enough density, and which are not in any. It has no k, it finds arbitrarily shaped clusters, and it explicitly labels points that belong to nothing as noise — scikit-learn returns those with the label -1. For zone design that noise label is a feature: those are the stops that need a separate long-run route rather than being forced into somebody’s round.

The two parameters

DBSCAN takes eps, a radius, and min_samples, a count. A point is a core point if at least min_samples points (including itself) lie within eps of it. Core points within eps of each other join the same cluster; non-core points within eps of a core point join as border points; everything else is noise. That is the entire algorithm, and both parameters do very different jobs.

min_samples controls how much evidence you demand before calling something a density region. The conventional floor for two-dimensional data is 4, and raising it makes the result more conservative and more resistant to a handful of stray points bridging two clusters. eps is the parameter people burn a day on, and it is the one that can be derived.

Deriving epsilon from stop density

For points scattered at random with density lambda per unit area — a homogeneous Poisson process — the expected distance from a point to its nearest neighbour is 1 / (2 × sqrt(lambda)). This is the Clark and Evans result from 1954, and it gives a principled anchor: eps should be a small multiple of that distance, because you want epsilon to comfortably reach a typical neighbour but not to reach across genuinely empty ground.

Take a concrete service area: 2,000 stops spread over 40 km².

lambda = 2000 / 40 km^2 = 50 stops per km^2 = 5.0e-5 stops per m^2

expected nearest-neighbour distance
  = 1 / (2 * sqrt(5.0e-5))
  = 1 / (2 * 0.0070711)
  = 1 / 0.0141421
  = 70.7 m

eps ~ 2x to 3x that  ->  140 m to 210 m; start at 175 m

Now sanity-check that against the physical situation rather than accepting it: 175 m is roughly two urban blocks, which is a plausible statement of “these stops belong to the same walking round”. If the number had come out at 8 m or at 4 km, the density input was wrong — usually because the bounding area included a large lake or an airport that contains no stops and inflates the denominator.

The complementary check is the k-distance plot: for every point, compute the distance to its min_samples-th nearest neighbour, sort those distances descending and plot them. The curve is flat and then bends sharply upward; the value at the knee is the epsilon that separates “typical neighbourhood” from “isolated”. When the derived figure and the knee agree you can stop tuning. When they disagree by a lot, the area is not uniformly dense, which is the real answer to the question and points you at HDBSCAN, whose whole purpose is clusters of differing density.

The metric trap

DBSCAN’s eps is in the units of whatever metric you pass it, and the default is Euclidean. Handing it raw latitude and longitude degrees means epsilon is in degrees and the space is anisotropic: at 51°N one degree of longitude is 62% of a degree of latitude, so your circular neighbourhood is an ellipse and clusters stretch east-west. The result looks superficially plausible, which is what makes it dangerous. This is the same class of bug as a distance calculation that comes out 20% wrong.

Two correct options. Either project to a local equal-distance coordinate reference system — a UTM zone, or a national grid — and cluster in metres, or use the haversine metric, in which case the inputs must be radians and epsilon must be an angle:

import numpy as np
from sklearn.cluster import DBSCAN

coords = np.radians(stops[["lat", "lon"]].to_numpy())
eps_rad = 175 / 6371008.8            # metres -> radians on the mean sphere

labels = DBSCAN(eps=eps_rad, min_samples=4,
                metric="haversine", algorithm="ball_tree").fit_predict(coords)

noise = (labels == -1).sum()

Forgetting the np.radians conversion does not raise: it silently treats degrees as radians, epsilon becomes meaningless, and everything collapses into one cluster. Check the noise count and the cluster count before you look at a map.

Density clusters are not zones

Here is where the obvious approach stops. DBSCAN puts no upper bound on cluster size, and in any real city one cluster will swallow the entire centre — three thousand stops that no single driver can serve. Density and workload are different quantities, and the algorithm only models the first.

• Treat DBSCAN as a first pass. Use it to separate the served area from the genuinely isolated stops, then subdivide any cluster whose total service time exceeds a shift into balanced sub-zones. Capacitated clustering, or a capacitated vehicle routing solver run directly on the cluster, is the tool for that step.
• Weight by service time, not by stop count. Twenty residential drops and twenty office-tower drops with a lift and a reception desk are not the same workload, so a balanced split on stop counts is unbalanced in practice.
• Cluster on network distance where it matters. Two stops 150 m apart across a motorway or a river are ten minutes apart. If barriers dominate your geography, a straight-line metric will merge zones that a driver cannot cross, and the fix is to cluster on a driving-time matrix instead.
• Zones must be stable. Re-running DBSCAN nightly on slightly different stops moves the boundaries, and drivers lose the route knowledge that makes them fast. Recompute on a slow cadence and assign new stops to existing zones in between.

What comes out of that is a set of zones; what happens inside one is a routing problem with its own complexity story.