Geostatistical Modeling Spatial Statistics with Python

Ripley's K Function Implementation Guide

TL;DR — Call ripley_k(coords, distances, window_polygon, n_simulations=999) with a NumPy (N, 2) projected coordinate array, a 1-D distance vector, and a Shapely Polygon observation window. The function returns k_obs (empirical K values) and k_sim (simulation matrix). Compare k_obs against np.pi * distances**2 to detect clustering (k_obs > πd²) or dispersion (k_obs < πd²).

Why This Matters

Ripley’s K function is the workhorse second-order statistic in point pattern analysis, evaluating whether point interactions — clustering, inhibition, or randomness — exist across a continuous range of spatial scales. Unlike the nearest-neighbour index, which collapses everything to a single scalar, K(d) exposes scale-dependent structure: a process may cluster at neighbourhood scale (d < 200 m) while exhibiting inhibition at micro-scale (d < 20 m) due to territorial constraints. This makes it indispensable for ecological species distribution studies, epidemiological hotspot detection, and urban infrastructure optimisation. The technique sits within the broader framework of Core Concepts of Spatial Statistics & Geostatistics as the canonical test for second-order spatial dependence, complementing first-order intensity surfaces from kernel density estimation and attribute-based tests like spatial autocorrelation metrics.

Environment and Version Pinning

python
# Install exact versions for reproducibility
# pip install numpy==1.26.4 scipy==1.13.1 shapely==2.0.4 matplotlib==3.9.0
python
import numpy as np          # 1.26.4
from scipy.spatial import cKDTree          # scipy 1.13.1
from shapely.geometry import Polygon, Point, MultiPoint  # shapely 2.0.4
import matplotlib.pyplot as plt            # 3.9.0

CRS requirement: All coordinates must be in a projected system (e.g., EPSG:32618 UTM Zone 18N, or a local planar CRS). The distances array must share the same linear unit as the coordinates (metres if UTM). Geographic coordinates (decimal degrees) will silently produce incorrect K values because the πd² baseline assumes Euclidean geometry.

Step-by-Step Implementation

Step 1 — Define the Observation Window

The observation window must tightly bound the area where events could theoretically occur. An over-large window inflates the denominator area |W| and suppresses the estimated intensity λ, biasing K upward (false clustering).

python
from shapely.geometry import Polygon

# Real use: load from a GeoDataFrame or GeoJSON
# window = gpd.read_file("study_area.gpkg").geometry.union_all()

# Synthetic example: 100 × 100 m square (projected, metres)
window = Polygon([(0, 0), (100, 0), (100, 100), (0, 100)])

# Quick sanity checks
assert window.is_valid, "Window polygon must be topologically valid"
assert window.area > 0, "Window must have positive area"

Step 2 — Compute Empirical K(d) with Isotropic Edge Correction

Raw neighbor counts suffer from boundary truncation: points near the window edge have fewer visible neighbors because the search radius extends outside the study area. The isotropic correction weights each focal point by the proportion of its search circle that lies inside the window.

python
def ripley_k(coords, distances, window_polygon, n_simulations=999):
    """
    Compute Ripley's K function with isotropic edge correction and
    a Monte Carlo CSR simulation envelope.

    Parameters
    ----------
    coords : np.ndarray, shape (N, 2)
        Projected point coordinates in metres.
    distances : np.ndarray, shape (M,)
        Evaluation distances (same unit as coords).
    window_polygon : shapely.geometry.Polygon
        Study area boundary in the same CRS as coords.
    n_simulations : int
        CSR simulations for the envelope. Use >=999 for publication.

    Returns
    -------
    k_obs : np.ndarray, shape (M,)   Observed K(d)
    k_sim : np.ndarray, shape (n_simulations, M)   Simulated K(d) matrix
    """
    n = len(coords)
    area = window_polygon.area
    lam = n / area          # global intensity: events per unit area
    tree = cKDTree(coords)

    # --- Observed K with isotropic (Ripley) edge correction ---
    k_obs = np.zeros(len(distances))

    for j, d in enumerate(distances):
        total = 0.0
        for i, pt in enumerate(coords):
            # Proportion of the search circle that falls inside the window
            circle = Point(pt).buffer(d, resolution=32)
            clipped_area = circle.intersection(window_polygon).area
            if clipped_area < 1e-12:
                continue
            # Edge weight: upscale count to compensate for hidden sector
            edge_weight = (np.pi * d ** 2) / clipped_area

            count = len(tree.query_ball_point(pt, d)) - 1  # exclude self
            total += count * edge_weight

        k_obs[j] = total / (lam ** 2 * area)

    # --- CSR simulation envelope (homogeneous Poisson, same n and window) ---
    xmin, ymin, xmax, ymax = window_polygon.bounds
    k_sim = np.zeros((n_simulations, len(distances)))
    rng = np.random.default_rng(42)

    for sim in range(n_simulations):
        # Rejection sampling: uniform candidates filtered to window interior
        sim_pts = []
        while len(sim_pts) < n:
            candidates = np.column_stack([
                rng.uniform(xmin, xmax, n),
                rng.uniform(ymin, ymax, n),
            ])
            inside = np.array([
                window_polygon.contains(Point(c)) for c in candidates
            ])
            sim_pts.extend(candidates[inside].tolist())
        sim_pts = np.array(sim_pts[:n])

        sim_tree = cKDTree(sim_pts)
        for j, d in enumerate(distances):
            counts = np.array([
                len(sim_tree.query_ball_point(p, d)) - 1 for p in sim_pts
            ])
            # No edge correction for simulated data (they already fill the window)
            k_sim[sim, j] = counts.sum() / (lam ** 2 * area)

    return k_obs, k_sim

Step 3 — Generate Clustered Test Data and Run

python
# Synthetic two-cluster pattern + random background
rng = np.random.default_rng(42)
pts = np.vstack([
    rng.normal([30, 30], 4, (25, 2)),   # cluster A
    rng.normal([70, 70], 4, (25, 2)),   # cluster B
    rng.uniform(0, 100, (20, 2)),       # background
])

# Clip to window to remove escaped points
xmin, ymin, xmax, ymax = window.bounds
mask = (
    (pts[:, 0] > xmin) & (pts[:, 0] < xmax) &
    (pts[:, 1] > ymin) & (pts[:, 1] < ymax)
)
pts = pts[mask]

distances = np.linspace(1, 30, 60)
k_obs, k_sim = ripley_k(pts, distances, window, n_simulations=199)

Step 4 — Apply the L-Function Transformation

The theoretical K baseline is K(d) = πd², a quadratic curve that amplifies visual variance at large distances. The L-function linearises the expectation and centres it at zero:

L(d)=K(d)π−dL(d) = \sqrt{\frac{K(d)}{\pi}} - d

Under complete spatial randomness, L(d)=0L(d) = 0. Values above zero indicate clustering; values below zero indicate regularity.

python
def l_transform(k_values, distances):
    """Apply the L-function variance-stabilising transform."""
    return np.sqrt(k_values / np.pi) - distances

l_obs = l_transform(k_obs, distances)
l_sim = l_transform(k_sim, distances)
l_lo  = l_sim.min(axis=0)
l_hi  = l_sim.max(axis=0)

Step 5 — Visualise the Results

python
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Panel 1: Raw K(d)
k_lo = k_sim.min(axis=0)
k_hi = k_sim.max(axis=0)
axes[0].fill_between(distances, k_lo, k_hi, alpha=0.2, color="gray",
                     label=f"CSR Envelope ({k_sim.shape[0]} sims)")
axes[0].plot(distances, k_obs, color="#e05c00", linewidth=2.5, label="Observed K(d)")
axes[0].plot(distances, np.pi * distances ** 2, "k--", label="Theoretical CSR (πd²)")
axes[0].set_xlabel("Distance d (m)")
axes[0].set_ylabel("K(d)")
axes[0].set_title("Ripley's K Function")
axes[0].legend(loc="upper left")
axes[0].grid(alpha=0.3)

# Panel 2: L-function
axes[1].fill_between(distances, l_lo, l_hi, alpha=0.2, color="gray",
                     label="CSR Envelope")
axes[1].plot(distances, l_obs, color="#e05c00", linewidth=2.5, label="Observed L(d)")
axes[1].axhline(0, color="black", linestyle="--", label="CSR baseline (L=0)")
axes[1].set_xlabel("Distance d (m)")
axes[1].set_ylabel("L(d)")
axes[1].set_title("L-Function (Variance-Stabilised)")
axes[1].legend(loc="upper left")
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig("ripley_k_output.png", dpi=150)
plt.show()

Visualising the K-Function Workflow

The diagram below illustrates how raw coordinates flow through edge correction, simulation, and the L-function transform to produce the final inference.

Ripley's K Function Computation Workflow Five stages: projected coordinates, isotropic edge correction, CSR Monte Carlo envelope, L-function transform, and spatial inference. Arrows connect each stage left to right. Projected Coordinates EPSG:32618 + Window Polygon Isotropic Edge Correction circle ∩ window weight = πd² / clipped area CSR Monte Carlo Envelope Poisson process ≥999 simulations min/max bands L-Function Transform √(K/π) − d baseline = 0 linear variance Spatial Inference L > 0 → clustering L < 0 → dispersion L ≈ 0 → random at scale d Step 1 Step 2 Step 3 Step 4 Step 5

Interpreting the Output

Output What it Means
k_obs[j] > np.pi * distances[j]**2 Clustering at scale distances[j]: more neighbors observed than random expectation
k_obs[j] < np.pi * distances[j]**2 Inhibition / regularity: inter-point repulsion or packing constraints
k_obs inside the CSR envelope No statistically detectable departure from complete spatial randomness
l_obs strongly positive at small d, flipping negative at large d Multi-scale structure — micro-clustering with macro-regularity (common in plant ecology)
Wide CSR envelope even with 999 sims Low point count (N < 50): K-function has high sampling variance; consider reporting the L-function instead

Critical Best Practices

Always Use a Projected CRS

Distances in the πd² formula are Euclidean. Using geographic coordinates (decimal degrees) makes the theoretical baseline wrong by orders of magnitude for anything beyond a few kilometres. Reproject with geopandas: gdf = gdf.to_crs("EPSG:32618") before extracting coordinates.

Match the Window to Actual Study Area

The global intensity λ = N / |W| is the denominator of K. If the window includes large uninhabitable regions (water bodies, restricted zones) where events cannot occur, λ is underestimated and K is inflated. Use the ecologically or administratively valid boundary, not a bounding box.

Use 999+ Simulations for Publication

With 99 simulations, the rank-based envelope test has an effective significance level of approximately 2/(99+1) = 0.02. For p < 0.005, use 999 simulations (alpha = 2/1000 = 0.002). The 199-simulation default in the example gives alpha ≈ 0.01.

Switch to the Inhomogeneous K-Function for Non-Stationary Processes

When the kernel density intensity surface shows strong gradients (e.g., higher event density near roads, rivers, or population centres), the homogeneous K-function will report spurious clustering even if no genuine inter-point interaction exists. The inhomogeneous K-function normalises by local intensity rather than global area. The PySAL pointpats library provides pointpats.Kest with inhomogeneous weighting. See also Stationarity & Trend Analysis for tests that help decide whether a homogeneous or inhomogeneous model is appropriate.

Scale Performance for Large Datasets

The per-point inner loop above runs in O(N × M) time where M is the number of distance bands. For N > 5,000:

  • Replace the per-simulation query_ball_point loop with tree.query_ball_tree(sim_tree, r=d) to get all pairs in one call.
  • For rectangular windows, replace the Shapely circle.intersection(window) with a closed-form rectangle–circle overlap formula (constant time per point, no geometry construction).
  • Parallelise across distance bands with joblib.Parallel(n_jobs=-1) — each band is independent.

Troubleshooting

Symptom Likely Cause Fix
k_obs everywhere above the CSR envelope, even at d=1 Window too large relative to the actual event domain Tighten window_polygon to the true observation area
k_obs below the theoretical πd² at all scales Geographic CRS used instead of projected Reproject to a planar CRS before running
CSR envelope extremely wide N < 30 — insufficient points for stable K estimation Consider quadrat methods or nearest-neighbour statistics instead
clipped_area consistently near zero for boundary points Window polygon uses a very different scale/units from coords Confirm both are in the same CRS and unit
Simulation runs very slowly (>10 min) Rejection sampling inefficient for complex window shapes Pre-compute a rasterised mask of the window and sample from valid pixels
L-function shows clustering at all scales with no plateau Inhomogeneous first-order intensity not accounted for Switch to the inhomogeneous K-function or detrend the intensity surface first

Next Steps

For a broader treatment of distance-based and quadrat statistics, return to the Point Pattern Analysis workflow page. If your study dataset has known spatial sampling bias — for example, opportunistic species records or mobile sensor pings — read Correcting Spatial Sampling Bias with GeoPandas before running K-function tests, since unweighted biased data will inflate clustering signals.

Related

← Back to Point Pattern Analysis