Point Pattern Analysis: Python Workflows for Spatial Point Data

Point pattern analysis (PPA) quantifies whether a set of discrete spatial events — disease cases, tree locations, crime incidents, sensor alerts — is randomly scattered, systematically clustered, or more regular than chance would predict. It belongs to the Core Concepts of Spatial Statistics & Geostatistics family of methods and is the right tool whenever your observations are exact coordinates rather than aggregate counts in polygons or values on a continuous raster surface.

Prerequisites #

• Python ≥ 3.10
• geopandas ≥ 0.14, pointpats ≥ 2.3, libpysal ≥ 4.9
• scipy ≥ 1.11, numpy ≥ 1.26, matplotlib ≥ 3.8
• Point data in a projected CRS (e.g., UTM, EPSG:326xx / 327xx, or a local state plane). Geographic lat/lon coordinates (WGS 84) will distort every distance metric and must be reprojected before analysis.
• An explicit study-area polygon defining where events could have been recorded. Without a defined observation window, edge-correction and intensity estimates are meaningless.
• Install: pip install geopandas pointpats libpysal scipy numpy matplotlib

Mathematical Core #

Point pattern analysis decomposes spatial structure into two levels that must be evaluated in order.

First-order intensity $\lambda(\mathbf{s})$ describes the expected number of events per unit area at location $\mathbf{s}$:

$$\lambda(\mathbf{s}) = \lim_{|d\mathbf{s}| \to 0} \frac{E[N(d\mathbf{s})]}{|d\mathbf{s}|}$$

When intensity is constant across the window ($\lambda(\mathbf{s}) = \lambda$), the process is homogeneous; otherwise it is inhomogeneous and driven by spatial covariates such as elevation, road proximity, or population density.

Second-order statistics characterise inter-point interaction after accounting for intensity. Ripley’s $K$ function is the canonical tool:

$$K(r) = \frac{1}{\lambda} \mathbb{E}\left[\text{number of additional events within distance } r \text{ of a randomly chosen event}\right]$$

Under Complete Spatial Randomness (CSR), $K(r) = \pi r^2$. The linearised form $L(r) = \sqrt{K(r)/\pi}$ centres the CSR baseline at $L(r) = r$, making deviations visually obvious:

$$L(r) - r > 0 \Rightarrow \text{clustering at scale } r$$

$$L(r) - r < 0 \Rightarrow \text{regularity (inhibition) at scale } r$$

The Nearest-Neighbour Index (NNI) provides a single-number global summary. Under CSR, the expected mean nearest-neighbour distance for $n$ points in area $A$ is:

$$\bar{d}_{\text{CSR}} = \frac{1}{2\sqrt{\lambda}}, \quad \lambda = \frac{n}{A}$$

$$\text{NNI} = \frac{\bar{d}_{\text{observed}}}{\bar{d}_{\text{CSR}}}$$

An NNI below 1 indicates clustering; above 1 indicates dispersion.

Step 1: Data Ingestion and Observation Window Definition #

A robust workflow begins with strict validation. Points must be in a metric projection and the observation window must be explicitly defined before any statistic is computed.

import geopandas as gpd
import numpy as np
from shapely.ops import unary_union

# Load point data and study boundary
points_gdf = gpd.read_file("event_locations.gpkg")
study_area  = gpd.read_file("study_boundary.gpkg")

# Enforce a projected CRS — UTM Zone 18N as an example
TARGET_CRS = "EPSG:32618"
points_gdf = points_gdf.to_crs(TARGET_CRS)
study_area  = study_area.to_crs(TARGET_CRS)

# Clip: keep only points that fall inside the window
window_geom = unary_union(study_area.geometry)
points_gdf  = points_gdf[points_gdf.geometry.within(window_geom)].copy()
points_gdf  = points_gdf.drop_duplicates(subset="geometry")  # remove exact duplicates

# Extract raw coordinate array (n × 2) for scipy / pointpats
coords = np.column_stack([points_gdf.geometry.x, points_gdf.geometry.y])
n      = len(coords)
area   = window_geom.area          # in square metres for EPSG:32618
lambda_hat = n / area              # global intensity: events per m²

print(f"Events: {n}  |  Area: {area/1e6:.2f} km²  |  λ̂: {lambda_hat*1e6:.4f} per km²")

The drop_duplicates call matters: coincident coordinates inflate nearest-neighbour counts and produce an artefactual NNI below 1 even in a genuinely random process.

Step 2: First-Order Intensity and Exploratory Mapping #

Before testing for inter-point interaction, characterise the baseline intensity. A kernel density estimate (KDE) reveals whether events are uniformly distributed or whether spatial gradients (population, infrastructure, ecology) concentrate them.

import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

# KDE with Scott's bandwidth — adjust bw_method for domain-specific scales
kde = gaussian_kde(coords.T, bw_method="scott")

# Evaluate on a 200×200 grid spanning the bounding box
xmin, ymin, xmax, ymax = points_gdf.total_bounds
gx, gy = np.mgrid[xmin:xmax:200j, ymin:ymax:200j]
intensity = kde(np.vstack([gx.ravel(), gy.ravel()])).reshape(200, 200)

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

study_area.plot(ax=axes[0], facecolor="none", edgecolor="black", lw=1.5)
points_gdf.plot(ax=axes[0], markersize=4, alpha=0.6, color="steelblue")
axes[0].set_title("Raw Point Distribution")

im = axes[1].imshow(
    intensity.T,
    extent=[xmin, xmax, ymin, ymax],
    origin="lower", cmap="viridis", aspect="auto"
)
study_area.plot(ax=axes[1], facecolor="none", edgecolor="white", lw=1)
plt.colorbar(im, ax=axes[1], label="Estimated intensity (events/m²)")
axes[1].set_title("KDE Intensity Surface")

plt.tight_layout()
plt.savefig("intensity_map.png", dpi=150, bbox_inches="tight")

If the KDE surface shows strong gradients, the process is inhomogeneous. Applying homogeneous second-order statistics to inhomogeneous data produces false positives for clustering — the pattern appears clustered simply because intensity is higher in one corner. In that case, proceed to Step 4 before drawing conclusions. For a detailed treatment of intensity gradients and trend removal, see the guide on stationarity and trend analysis.

Step 3: Nearest-Neighbour Index #

The NNI offers a quick global diagnostic before committing to the more computationally intensive K-function analysis.

from scipy.spatial import cKDTree

tree = cKDTree(coords)

# k=2 because k=1 returns the point itself (distance = 0)
nn_dists, _ = tree.query(coords, k=2)
nn_dists     = nn_dists[:, 1]   # first *other* neighbour

expected_nn  = 1.0 / (2.0 * np.sqrt(lambda_hat))
observed_nn  = nn_dists.mean()
nni          = observed_nn / expected_nn

# Z-score for large samples (n > 30)
se_csr = np.sqrt((4.0 - np.pi) / (4.0 * np.pi * n * lambda_hat))
z      = (observed_nn - expected_nn) / se_csr

print(f"Observed mean NN : {observed_nn:.2f} m")
print(f"Expected mean NN : {expected_nn:.2f} m (CSR)")
print(f"NNI              : {nni:.4f}  (<1 clustering, >1 regularity)")
print(f"Z-score          : {z:.3f}   (|z| > 1.96 → p < 0.05)")

Interpret carefully: the NNI is a global average and will return a value near 1 for a pattern that has both tight sub-clusters and large empty areas — the clustering and the voids cancel out. Use it as a screening tool, not a final answer.

Step 4: Scale-Dependent Clustering with Ripley’s K and L #

Ripley’s K reveals at which spatial scales clustering or regularity operates, which the NNI cannot. The pointpats library provides edge-corrected estimators.

from pointpats import PointPattern, K_test, ripley_k

# pointpats expects coordinates as an (n, 2) array
pp = PointPattern(coords)

# Define distance bands — typically up to half the shortest study-area dimension
r_max   = min((xmax - xmin), (ymax - ymin)) / 4   # conservative upper bound
r_vals  = np.linspace(0, r_max, 50)

# Ripley's K with Ripley's edge correction ("r" kernel)
k_obs = ripley_k(pp, support=r_vals, keep_simulations=False)

# Transform to L(r) and centre: L(r) - r
l_obs = np.sqrt(k_obs / np.pi) - r_vals

plt.figure(figsize=(8, 4))
plt.axhline(0, color="grey", linestyle="--", label="CSR expectation")
plt.plot(r_vals, l_obs, color="#7c3aed", lw=2, label="Observed L(r) − r")
plt.xlabel("Distance r (m)")
plt.ylabel("L(r) − r")
plt.title("Ripley's L function — deviation from CSR")
plt.legend()
plt.tight_layout()
plt.savefig("ripley_l.png", dpi=150, bbox_inches="tight")

For publication-quality inference, add Monte Carlo simulation envelopes — the full guide with envelope construction and significance bands is in the Ripley’s K Function Implementation Guide.

Diagnostic Configuration: Bandwidth and Distance Band Selection #

Two parameters drive the most common errors in PPA.

KDE bandwidth (bw_method in gaussian_kde):

• Scott’s rule ("scott") is safe for a first pass but assumes a roughly Gaussian density.
• Use likelihood cross-validation (pointpats.kde.LikelihoodCVBandwidth) when the intensity surface will feed into inhomogeneous K-function corrections.
• Domain knowledge should constrain the plausible range: a bandwidth of 500 m for a city-block crime dataset is far too coarse; 50 m may be right.

Ripley’s K distance range (r_max):

• The standard rule-of-thumb caps $r$ at one-quarter of the shorter study-area dimension to keep edge effects manageable.
• Never extend $r$ beyond half the shortest dimension; beyond that, edge correction assumptions break down and variance explodes.
• For ecological data with known dispersal distances (e.g., seed dispersal range), set r_max to that domain-specific upper bound.

# Example: likelihood cross-validation bandwidth via pointpats
from pointpats.kde import LikelihoodCVBandwidth

bw_selector = LikelihoodCVBandwidth(coords)
bw_opt      = bw_selector.bandwidth          # optimal h in map units
print(f"Cross-validated bandwidth: {bw_opt:.1f} m")

kde_cv = gaussian_kde(coords.T, bw_method=bw_opt / coords.std())

Step 5: Handling Inhomogeneous Processes #

Most real-world point processes are inhomogeneous — event intensity varies because of underlying drivers (population density, habitat suitability, network proximity) rather than genuine clustering. Applying homogeneous K-function tests to such data confounds environmental gradients with true interaction. The corrective approach uses inverse-intensity weights when building the spatial weight matrices that feed residual analysis, or feeds directly into an inhomogeneous K estimator.

# Step 1: Estimate intensity at each observed point (leave-one-out KDE)
kde_global        = gaussian_kde(coords.T, bw_method="scott")
intensity_at_pts  = kde_global(coords.T)

# Step 2: Inverse-intensity weights — down-weight high-intensity regions
weights     = 1.0 / intensity_at_pts
weights    /= weights.mean()          # normalise so mean weight = 1

# Step 3: Effective sample size (Kish approximation) — confirms weighting is sane
ess = n / ((weights**2).sum() / n)
print(f"Effective sample size after IIW: {ess:.1f}  (raw n={n})")

# Weights >> 1 mark points in sparse regions that contribute more to
# residual analysis; weights << 1 mark points in dense 'expected' areas.

The inhomogeneous K estimator in pointpats accepts these weights directly via the weights argument of ripley_k, producing a statistic that tests for interaction beyond the trend. This is the correct approach when environmental covariates are available — analogous to using sampling bias mitigation techniques that correct for unequal detection probabilities across space.

Output Interpretation #

The L function curve should be interpreted alongside its Monte Carlo envelope, not in isolation. A curve that strays outside the envelope at only one or two distance bands may be a Type I error — especially if you evaluate many $r$ values without a multiple-comparisons correction (use a global rank envelope test from spatstat-style methods).

Production Considerations #

Scaling to large point sets. Exact pairwise distance matrices are $O(n^2)$ in memory and time. For $n > 50{,}000$:

• Use scipy.spatial.cKDTree for $O(n \log n)$ nearest-neighbour queries.
• Switch to grid-based quadrat counting (divides the window into cells and applies a chi-square test) for a fast first-order check.
• For K-function estimation, use approximate methods in pointpats or subsample the point set to $n \approx 5{,}000$ and rescale $\lambda$.

Reproducibility. Monte Carlo envelopes require a fixed random seed: np.random.default_rng(42). Store the seed, the window geometry (as WKT), and the library versions alongside your results.

Multi-type patterns. When events belong to two or more types (e.g., predator vs. prey locations), use the cross K-function $K_{12}(r)$ which tests whether type-1 events are clustered around type-2 events beyond chance. pointpats provides cross_k.

Memory limits. For windows with complex boundaries (high-vertex polygons), edge correction at each distance band requires point-in-polygon tests that dominate runtime. Simplify the boundary with window_geom.simplify(tolerance) before passing to pointpats.

For patterns with strong covariates, extend to a Poisson regression intensity model (log-linear in covariates) and use the fitted intensity surface as the weights input — this ties directly into the spatial autocorrelation metrics workflow where residual Moran’s I can confirm whether residual interaction remains after trend removal.

Troubleshooting #

Next Steps #

Once you can confidently distinguish clustering from inhomogeneity, the logical next step is to quantify clustering across a range of scales using Monte Carlo simulation envelopes — the Ripley’s K Function Implementation Guide walks through the full edge-corrected implementation with significance bands. For workflows where spatial bias in data collection (citizen science, mobile sensors) may have generated the apparent pattern, see Correcting Spatial Sampling Bias with GeoPandas.

Frequently Asked Questions #

What is Complete Spatial Randomness and why does it matter?

CSR (homogeneous Poisson process) is the null model for point pattern analysis. Under CSR every location in the window is equally likely to receive an event and points do not attract or repel each other. Every clustering or dispersion test — NNI, Ripley’s K, quadrat chi-square — measures deviation from this baseline. Failing to account for inhomogeneous intensity inflates apparent clustering because the CSR assumption is violated.

What CRS should I use for point pattern analysis?

Always reproject to a local projected CRS (UTM, state plane) before computing distances or areas. Geographic coordinates (WGS 84 lat/lon) produce angular distances that are non-uniform in metres and will distort every distance-based statistic.

How do I choose a bandwidth for KDE intensity estimation?

Scott’s rule works as a starting point but for spatial intensity you typically need cross-validated or domain-informed bandwidths. A bandwidth too small produces a spiky surface that masks regional patterns; too large over-smooths and hides local concentrations. Use likelihood cross-validation from pointpats or run KDE at several bandwidths and compare with quadrat variance tests.

When should I use Ripley’s K versus the nearest-neighbour index?

The NNI gives a single global summary; it misses multi-scale structure. Ripley’s K evaluates clustering at every distance band simultaneously, so it can reveal that a pattern is clustered at short ranges but regular at long ranges. Use the NNI for a quick diagnostic and Ripley’s K for any decision that requires publication-quality evidence.

Related #

• Ripley’s K Function Implementation Guide — edge-corrected estimators, Monte Carlo envelopes, and significance testing in Python
• Stationarity and Trend Analysis — detecting and removing first-order trends before second-order testing
• Spatial Autocorrelation Metrics — Moran’s I and LISA for areal and lattice data
• Sampling Bias Mitigation — correcting for unequal detection probabilities that mimic clustering

← Back to Core Concepts of Spatial Statistics & Geostatistics