Nearest Neighbour vs Ripley's K for Clustering

TL;DR: Use the Clark-Evans index — pointpats.PointPattern(coords).mean_nnd divided by 0.5 / np.sqrt(n / area) — when one characteristic scale is enough. Use pointpats.l_test(coords, support=..., hull=window, linearized=True, n_simulations=199) when the scale of clustering is the question. R reads each point’s nearest neighbour only; L reads every distance band.

Why This Matters #

Both tests answer “is this pattern clustered?” against the same null hypothesis of complete spatial randomness, and on a well-behaved single-scale process they agree. The difference is what each throws away. The nearest-neighbour index keeps one distance per point — the shortest — and averages it, so the pattern is described at whatever scale that mean happens to sit at. Ripley’s K keeps every inter-point distance and reports the expected number of further points within radius $d$ of a typical point, for every $d$ you ask for. One is a number; the other is a curve, and a curve is the only thing that can say at what distance the process departs from randomness.

That distinction decides the shape of the rest of the analysis. If clustering turns out to be confined to 200–400 metres, that number sets the bandwidth for a density surface and the lag spacing for a variogram; a single R of 0.69 sets nothing. This page sits inside Point Pattern Analysis and leans on the same distributional machinery as the rest of Core Concepts of Spatial Statistics & Geostatistics.

The two statistics are built from the same intensity $\lambda = n / |A|$, where $|A|$ is the area of the observation window. Under complete spatial randomness the mean distance from a point to its nearest neighbour is

$$\mathbb{E}[\bar{d}] = \frac{1}{2\sqrt{\lambda}}, \qquad \mathrm{SE}[\bar{d}] = \sqrt{\frac{4 - \pi}{4\pi n \lambda}} \approx \frac{0.26136}{\sqrt{n\lambda}},$$

and the Clark-Evans ratio is $R = \bar{d}_{\text{obs}} / \mathbb{E}[\bar{d}]$, with $R < 1$ indicating clustering, $R > 1$ regularity, and $R = 1$ the null. Ripley’s K instead counts:

$$\hat{K}(d) = \frac{1}{\lambda n} \sum_{i} \sum_{j \neq i} \mathbf{1}\left( \lVert s_i - s_j \rVert \le d \right), \qquad \mathbb{E}[K(d)] = \pi d^2 .$$

The variance-stabilising transform $L(d) = \sqrt{K(d)/\pi}$ turns the parabola $\pi d^2$ into the straight line $L(d) = d$, so plotting $L(d) - d$ puts the null at zero and makes departures readable at every distance on the same vertical scale. That is the form worth plotting; the raw $K$ curve is dominated by its own quadratic growth.

Environment and Version Pinning #

pointpats supplies both the nearest-neighbour distances and the distance-function tests. shapely carries the observation window, and scipy is only needed for the normal tail probability behind the Clark-Evans z-score.

pip install "pointpats==2.5.0" "libpysal==4.12.1" "geopandas==1.0.1" \
            "shapely==2.0.4" "numpy==1.26.4" "scipy==1.13.1"

import numpy as np                          # 1.26.4
import geopandas as gpd                     # 1.0.1
from shapely.geometry import Polygon        # 2.0.4
from scipy.stats import norm                # 1.13.1
from pointpats import PointPattern, l_test  # 2.5.0

Step-by-Step Implementation #

1. Fix the window before anything else #

Both statistics divide by the intensity $\lambda = n / |A|$, so the answer is a function of the polygon you call the study area. Take it from the sampling design — the surveyed forest compartment, the licensed quarry block, the borough boundary — not from the convex hull of the points, which shrinks towards the data and inflates $\lambda$.

SIDE = 2000.0                                     # metres
WINDOW = Polygon([(0, 0), (SIDE, 0), (SIDE, SIDE), (0, SIDE)])
AREA = WINDOW.area                                # 4_000_000 m^2

# In production the window and the points arrive as layers:
#   pts = gpd.read_file("events.gpkg").to_crs("EPSG:32630")
#   window = gpd.read_file("survey_area.gpkg").to_crs(pts.crs).union_all()
#   coords = np.column_stack([pts.geometry.x, pts.geometry.y])

2. Generate a pattern with two distinct scales #

The pattern below is deliberately awkward: offspring are scattered in discs of 200 m radius around 16 parents whose locations follow an environmental gradient rising to the north-east, and no point may fall within 45 m of another. That is short-range inhibition inside long-range clustering.

def two_scale_pattern(seed=7, n_parents=16, per_parent=25,
                      cluster_radius=200.0, hard_core=45.0, side=SIDE):
    """Gradient-driven patches with a hard-core exclusion inside them."""
    rng = np.random.default_rng(seed)
    cand = rng.uniform(0, side, size=(4000, 2))
    weight = ((cand[:, 0] + cand[:, 1]) / (2 * side)) ** 3
    parents = cand[rng.random(4000) < weight][:n_parents]

    pts = []
    for a, b in parents:
        placed, tries = 0, 0
        while placed < per_parent and tries < 20_000:
            tries += 1
            theta = rng.uniform(0, 2 * np.pi)
            r = cluster_radius * np.sqrt(rng.random())
            x, y = a + r * np.cos(theta), b + r * np.sin(theta)
            if not (0 <= x <= side and 0 <= y <= side):
                continue
            if pts:
                arr = np.asarray(pts)
                if np.hypot(arr[:, 0] - x, arr[:, 1] - y).min() < hard_core:
                    continue
            pts.append((x, y))
            placed += 1
    return np.asarray(pts)


coords = two_scale_pattern()
print(coords.shape)

(384, 2)

3. Compute the Clark-Evans index #

pointpats gives you the mean nearest-neighbour distance; it does not ship a Clark-Evans function, and that is just as well, because the ratio depends entirely on which area you divide by. Assemble it explicitly against the window from step 1.

def clark_evans(coords, area):
    """Clark-Evans R with its normal-approximation z-score and p-value."""
    pp = PointPattern(coords)
    n = len(coords)
    lam = n / area
    d_obs = pp.mean_nnd
    d_exp = 0.5 / np.sqrt(lam)
    se = 0.26136 / np.sqrt(n * lam)
    z = (d_obs - d_exp) / se
    return d_obs, d_exp, d_obs / d_exp, z, 2 * (1 - norm.cdf(abs(z)))


d_obs, d_exp, R, z, p = clark_evans(coords, AREA)
print(f"n = {len(coords)}, intensity = {len(coords) / AREA:.2e} points/m2")
print(f"observed mean nnd = {d_obs:.2f} m")
print(f"expected mean nnd = {d_exp:.2f} m")
print(f"Clark-Evans R = {R:.3f}   z = {z:.2f}   p = {p:.3f}")

n = 384, intensity = 9.60e-05 points/m2
observed mean nnd = 51.67 m
expected mean nnd = 51.03 m
Clark-Evans R = 1.012   z = 0.47   p = 0.639

The verdict is unambiguous and wrong: R of 1.012 with p = 0.639 says the pattern is indistinguishable from complete spatial randomness.

4. Run the L function with a simulation envelope #

l_test computes the statistic at each distance in support, then repeats it on n_simulations random patterns drawn inside hull. Passing linearized=True returns $L(d) - d$, and keep_simulations=True retains the simulated curves so you can build the envelope yourself.

support = np.arange(25.0, 801.0, 25.0)
res = l_test(coords, support=support, hull=WINDOW, linearized=True,
             n_simulations=199, keep_simulations=True)

lo = res.simulations.min(axis=0)
hi = res.simulations.max(axis=0)

for d in (25, 50, 100, 250, 500, 800):
    i = int(np.argmin(np.abs(res.support - d)))
    flag = ("below" if res.statistic[i] < lo[i]
            else "above" if res.statistic[i] > hi[i] else "inside")
    print(f"d={d:4.0f} m  L-d={res.statistic[i]:7.2f}  "
          f"envelope=[{lo[i]:7.2f},{hi[i]:7.2f}]  {flag}")

A representative run:

d=  25 m  L-d= -25.00  envelope=[  -4.64,   5.82]  below
d=  50 m  L-d=   5.44  envelope=[  -7.01,   6.06]  inside
d= 100 m  L-d=  48.56  envelope=[  -8.58,   2.47]  above
d= 250 m  L-d=  88.32  envelope=[ -19.97,  -7.44]  above
d= 500 m  L-d=  41.30  envelope=[ -66.52, -37.08]  above
d= 800 m  L-d= -60.88  envelope=[-161.22,-115.72]  above

At 25 m the observed value is exactly $-25$, because the hard core forbids every pair closer than 45 m and $K(25) = 0$. From 100 m to 800 m the curve is above the envelope, peaking at $+88.3$ at 250 m — roughly the diameter of a patch. The same 384 points are regular at short range and strongly clustered at medium range, and the one-number test saw neither.

5. Confirm that the two agree when there is only one scale #

None of this makes R useless. On a plain Thomas process — patches with no inhibition inside them — the mean nearest-neighbour distance is genuinely shortened, and both tests say clustered.

def thomas_pattern(seed=11, n_parents=20, per_parent=20, sigma=90.0, side=SIDE):
    rng = np.random.default_rng(seed)
    px, py = rng.uniform(0, side, n_parents), rng.uniform(0, side, n_parents)
    pts = []
    for a, b in zip(px, py):
        for _ in range(per_parent):
            while True:
                x, y = a + rng.normal(0, sigma), b + rng.normal(0, sigma)
                if 0 <= x <= side and 0 <= y <= side:
                    pts.append((x, y))
                    break
    return np.asarray(pts)


clustered = thomas_pattern()
d_obs, d_exp, R, z, p = clark_evans(clustered, AREA)
print(f"n = {len(clustered)}  observed = {d_obs:.2f} m  expected = {d_exp:.2f} m")
print(f"Clark-Evans R = {R:.3f}   z = {z:.2f}   p = {p:.3f}")

n = 400  observed = 34.34 m  expected = 50.00 m
Clark-Evans R = 0.687   z = -11.98   p = 0.000

R of 0.687 is a decisive rejection, obtained in milliseconds without a single simulation. Run l_test on the same pattern and the curve is above the envelope from 25 m to about 600 m — the extra information is the range, not the verdict.

Interpreting the Output #

Read R as a ratio and a z-score together. Values below about 0.8 with a large negative z mean the shortest links are much shorter than chance; values above about 1.2 with a large positive z mean the opposite. The trap is the middle: R near 1 is not evidence of randomness, it is the absence of evidence at one scale, and the two-scale pattern above is exactly how that arises.

Read the L curve as a sequence of excursions rather than a single verdict. Where it leaves the envelope upwards, points have more neighbours within that radius than randomness allows; where it leaves downwards, they have fewer. The distance at which the excursion begins is the inner scale of the process, the distance at which it peaks is roughly the patch diameter, and the distance at which it returns is the outer scale. Those are the three numbers worth quoting.

The envelope drifts below zero as $d$ grows, and that is expected, not a bug. Without edge correction the pair counts near the boundary are too low, so both the observed curve and every simulated curve are biased downward by the same amount. Comparing the observed curve with the theoretical zero line instead of with the envelope makes every pattern look regular at large $d$. This is also why l_test must receive the real window through hull: the simulations are drawn inside it, and an over-tight hull produces an over-tight envelope. The full estimator, its edge-correction options and its computational cost are worked through in the Ripley’s K-Function Implementation Guide.

Critical Best Practices #

Project first, and check the units are metres #

Both statistics are ratios involving points per unit area, so running either on longitude and latitude produces a number with no interpretation: degrees of longitude shrink with latitude, the window area is wrong, and R will drift with the study area’s position on the globe. Reproject to a UTM zone or a national grid before the first call, following reprojecting CRS for accurate distance calculations, and assert pts.crs.is_projected in the pipeline rather than trusting the file.

Never let the convex hull stand in for the window #

The default hull in pointpats is the bounding box of the points, and hull="convex" is tighter still. Both shrink towards the observed data, which raises $\lambda$, lowers the expected nearest-neighbour distance and pushes R upward — a clustered pattern can be dragged towards 1 by nothing more than an unstated default. Pass the surveyed polygon explicitly, and if the survey area genuinely is unknown, say so and report how much R moves between the bounding box and the convex hull.

Do not expect edge_correction to do anything #

The argument exists in the pointpats signature but is not implemented in the current release; supplying anything other than None raises NotImplementedError. This is less serious than it sounds, because the simulations are drawn in the same window and inherit the same downward bias, so the comparison remains valid. It does mean the observed curve is not an unbiased estimate of $K$, and that you cannot compare a curve computed on one window against a curve computed on another.

A pointwise envelope is not a global test #

Taking the minimum and maximum of 199 simulations gives a pointwise level of about 1% at each distance, but the curve is inspected at all 32 distances in the support simultaneously, so the chance of some excursion under a true null is much higher than 1%. Use the envelope to see the shape and read the scale; if a defensible p-value is needed, fix the distance range in advance and summarise the deviation over it, or use a rank envelope test.

Homogeneous K assumes a constant intensity #

The gradient in the example above is a first-order effect — more points in the north-east — while a hard core is a second-order effect between points. Ripley’s K cannot separate them, and reports both as clustering. If the intensity plainly varies, estimate it first, ideally with a surface built as described in kernel density estimation bandwidth selection in Python, and then move to the inhomogeneous K function, which weights each pair by the reciprocal of the local intensity.

Troubleshooting #

Next Steps #

Once the L curve has given you a scale, use it: it sets the plausible bandwidth range for a density surface and the maximum lag worth fitting in a variogram. Work through the estimator itself in the Ripley’s K-Function Implementation Guide, then return to Point Pattern Analysis for the wider set of second-order methods.

Frequently Asked Questions #

Can the nearest-neighbour index miss clustering that Ripley’s K detects?

Yes, and it happens whenever the process has more than one scale. The Clark-Evans index uses only each point’s single shortest link, so it describes the pattern at roughly the mean nearest-neighbour distance and nowhere else. A pattern with hard-core inhibition at 45 metres inside patches that are themselves clustered at 250 metres can return R of 1.01 with a p-value of 0.64, because the short-range regularity and the long-range clustering pull the one number in opposite directions and cancel.

Why does the simulation envelope sit below zero at long distances?

Because the estimator has no edge correction. Points near the window boundary have part of their neighbourhood outside the study area, so the pair counts are too low, and the bias grows with distance. Simulating complete spatial randomness inside the same window reproduces exactly the same bias, which is why the envelope drifts downward. Compare the observed curve with the envelope, never with the theoretical value of zero, or every pattern will look regular at large distances.

How many simulations does the envelope need?

For a pointwise envelope drawn as the minimum and maximum of the simulations, 199 runs give a nominal two-sided pointwise level of one per cent. That is enough to see the shape of the departure. It is not a valid global test, because the curve is inspected at every distance in the support at once. If you need a p-value you can defend, use a rank envelope or a maximum-absolute-deviation summary over a fixed distance range decided before looking at the data.

Does either test require a projected coordinate reference system?

Both do. The nearest-neighbour index divides an observed distance by an expectation derived from points per unit area, and Ripley’s K divides a pair count by the same intensity, so degrees of latitude and longitude produce a ratio with no meaning. Reproject to a metric CRS such as a UTM zone or a national grid before either call, and make sure the polygon used as the observation window carries the same CRS as the points.

Related

• Ripley’s K-Function Implementation Guide — the estimator, its edge corrections and its cost at scale
• Kernel Density Estimation Bandwidth Selection in Python — estimating the intensity surface the homogeneous K assumes away
• Reprojecting CRS for Accurate Distance Calculations — the metric coordinates both tests require

← Back to Point Pattern Analysis