Fixing Edge Effects in Poisson Point Patterns

Clustering statistics computed on a synthetic Poisson point pattern read as artificially clustered near the boundary, not because the process is clustered but because points close to the extent have no neighbours counted beyond it.

Part of Point Process Simulation Models: that workflow covers generating spatial point patterns; this page isolates the boundary-bias failure that corrupts second-order summary statistics — nearest-neighbour distances and Ripley’s K — near the edge of the study window, and the three standard remedies that restore an unbiased estimate.

Root Cause: The Observation Window Truncates Every Neighbourhood

A point process is defined on the whole plane, but you only ever observe it inside a finite window WW. Any statistic that counts neighbours within a distance rr of a point implicitly assumes the full disc of radius rr around that point is observed. For a point within rr of the boundary, part of that disc lies outside WW and its points are simply absent from the count. The estimate is biased downward — fewer neighbours than the process actually produced — and because the deficit grows as points approach the edge, the bias masquerades as a spatial gradient in intensity or clustering.

Ripley’s K function makes the mechanism explicit. For a stationary process of intensity λ\lambda, λK(r)\lambda K(r) is the expected number of further points within distance rr of a typical point. The naive estimator is

K^(r)=Wn(n1)iji1 ⁣(dijr)\hat{K}(r) = \frac{|W|}{n(n-1)} \sum_{i} \sum_{j \ne i} \mathbf{1}\!\left(d_{ij} \le r\right)

where W|W| is the window area and dijd_{ij} the inter-point distance. The indicator only ever counts pairs inside WW, so points near the boundary contribute too few pairs and K^(r)\hat{K}(r) falls below the theoretical Poisson value K(r)=πr2K(r) = \pi r^2. The relative bias grows with rr: once rr approaches the window’s short dimension, a large fraction of every disc spills outside and the estimate collapses.

The same truncation biases the nearest-neighbour distance distribution G(r)G(r): a point near the corner may have its true nearest neighbour just outside WW, so the observed nearest-neighbour distance is overstated, pushing G(r)G(r) toward apparent regularity. Both statistics are second-order — they depend on pairs, not marginal counts — which is exactly why the realism metrics used to score synthetic against source will flag an edge-biased pattern as failing to match even a correctly generated source.

Boundary bias in a Poisson point pattern and its correction by a guard region and toroidal wrapping Three panels. Left, labelled naive window: a square window with scattered points; an interior point carries a full circular search disc fully inside the window, while a point near the bottom-right corner carries a disc that is clipped by the two window edges, so neighbours outside are not counted and the count is biased low. Middle, labelled guard region: an outer generated buffer square contains an inner dashed analysis window inset by the search radius; a point on the inner window still has its full disc covered because the buffer supplies the neighbours. Right, labelled toroidal wrap: a point near the right edge has a search disc that exits the right edge and re-enters from the left edge, shown as two arcs, so the disc is always complete. Naive window full disc clipped → undercount Guard region generated buffer analysis window (inset by r) Toroidal wrap disc re-enters opposite edge
The naive count undercounts near the boundary; a guard region analyses only points whose full disc is observed, while toroidal wrapping supplies the missing neighbours from the opposite edge.

Minimal Reproducer: Show the Downward Bias

Generate a homogeneous Poisson pattern with a known intensity and compare the naive K^(r)\hat K(r) against the theoretical πr2\pi r^2. Pin the toolchain so distances and RNG streams match across machines.

numpy==1.26.*
scipy==1.11.*        # cKDTree for pair counting
python
import numpy as np
from scipy.spatial import cKDTree

def poisson_pattern(intensity: float, size: float, seed: int) -> np.ndarray:
    """Homogeneous Poisson process on a [0, size]^2 window."""
    rng = np.random.default_rng(seed)
    n = rng.poisson(intensity * size**2)
    return rng.uniform(0.0, size, size=(n, 2))

def naive_k(pts: np.ndarray, r: float, size: float) -> float:
    tree = cKDTree(pts)
    pairs = tree.count_neighbors(tree, r) - len(pts)   # subtract self-pairs
    lam = len(pts) / size**2
    return pairs / (lam * len(pts))                    # |W| n_hat cancels into lambda form

pts = poisson_pattern(intensity=500 / 1.0, size=1.0, seed=20260627)
for r in (0.05, 0.15, 0.25):
    print(r, round(naive_k(pts, r, 1.0), 4), "vs theoretical", round(np.pi * r**2, 4))
# naive K falls further below pi r^2 as r grows -> boundary bias, not clustering

The gap between naive_k and πr2\pi r^2 widens with rr: that widening deficit is the edge effect, and it is systematic, not sampling noise.

Fix: Guard Region, Toroidal Wrap, and an Edge-Corrected Estimator

Three remedies address the bias at different points in the pipeline. A guard region generates points in a buffer wider than the analysis window and only analyses points whose full search disc is observed. Toroidal wrapping treats the window as a torus so a disc that exits one edge re-enters the opposite one — valid only for a stationary process on a rectangular window. An edge-corrected estimator (Ripley’s isotropic correction) reweights each pair by the reciprocal of the observed fraction of its disc, correcting the statistic without changing the pattern.

python
import numpy as np
from scipy.spatial import cKDTree

def generate_with_guard(intensity: float, window: float, r_max: float, seed: int):
    """Generate on a buffered window; return (all_points, analysis_mask).
    Analysis points are those at least r_max from every boundary, so their
    full disc is covered by generated neighbours in the guard band."""
    buffered = window + 2 * r_max
    rng = np.random.default_rng(seed)
    n = rng.poisson(intensity * buffered**2)
    pts = rng.uniform(0.0, buffered, size=(n, 2)) - r_max   # shift so window is [0, window]
    inside = np.all((pts >= 0.0) & (pts <= window), axis=1)
    return pts, inside                                      # analyse pts[inside] against ALL pts

def toroidal_k(pts: np.ndarray, r: float, size: float) -> float:
    """Ripley's K under toroidal (periodic) boundary — no edge deficit by construction."""
    diff = np.abs(pts[:, None, :] - pts[None, :, :])
    diff = np.minimum(diff, size - diff)                    # wrap each axis onto the torus
    d = np.sqrt((diff**2).sum(axis=-1))
    np.fill_diagonal(d, np.inf)                             # drop self-pairs
    lam = len(pts) / size**2
    return (d <= r).sum() / (lam * len(pts))

For an irregular (non-rectangular) window where wrapping is invalid, use the isotropic edge correction: weight each ordered pair (i,j)(i,j) with dijrd_{ij} \le r by wijw_{ij}, the reciprocal of the proportion of the circle of radius dijd_{ij} centred at ii that lies inside WW.

python
import numpy as np

def isotropic_k(pts, r, window_area, edge_weight):
    """Ripley isotropic-corrected K. edge_weight(i, j) returns 1 / (observed disc fraction).
    Reweighting recovers the pairs the boundary hid, without moving any point."""
    n = len(pts)
    lam = n / window_area
    total = 0.0
    d = np.sqrt(((pts[:, None, :] - pts[None, :, :]) ** 2).sum(-1))
    for i in range(n):
        for j in range(n):
            if i != j and d[i, j] <= r:
                total += edge_weight(i, j)                  # up-weight boundary pairs
    return total / (lam * n)

Prefer the guard region when you control generation and can afford to synthesize a buffer — it is the cleanest because the statistic then needs no correction at all. Reach for the isotropic estimator when the window is fixed and irregular. The guard-region approach also composes naturally with a downstream polygon tessellation step, since the buffer points seed boundary cells that would otherwise be malformed.

Verification Step: Assert the Corrected K Matches Poisson

Under complete spatial randomness the corrected K^(r)\hat K(r) should track πr2\pi r^2 within a Monte Carlo envelope. Gate on the relative error at large rr, where the naive estimator fails hardest, and confirm the naive estimator still fails so the test cannot silently pass on a broken correction.

python
import numpy as np

def test_edge_correction_removes_bias():
    size, intensity, r = 1.0, 500.0, 0.25
    theo = np.pi * r**2

    corrected = np.mean([
        toroidal_k(poisson_pattern(intensity, size, s), r, size)
        for s in range(30)                                  # Monte Carlo over seeds
    ])
    # 1. Corrected estimator matches theory within 5% at large r.
    assert abs(corrected - theo) / theo < 0.05, f"corrected K biased: {corrected} vs {theo}"

    # 2. Naive estimator MUST still be biased low -> proves the test exercises the effect.
    naive = np.mean([
        naive_k(poisson_pattern(intensity, size, s), r, size)
        for s in range(30)
    ])
    assert naive < theo * 0.9, "naive K not biased -> reproducer no longer exercises edge effect"

The second assertion is the guard against a false pass: if a refactor accidentally applied the correction inside naive_k, the first assertion alone would still pass while the demonstration of the effect quietly disappeared.

Edge Cases & Gotchas

Toroidal wrapping on a non-stationary or irregular window. Wrapping assumes the process is stationary and the window tiles the plane. Apply it to a coastline-clipped study area or a pattern with an intensity gradient and you fabricate neighbour relationships across a boundary that has no physical meaning, biasing K upward. Restrict toroidal correction to rectangular windows and verified stationarity; otherwise use the guard region or isotropic weights.

Guard-band width smaller than r_max. If the buffer is narrower than the largest analysis radius, corner analysis points still have clipped discs and the bias returns for large rr only — a subtle failure that looks like clustering at scale. Size the guard band to at least r_max on every side, and assert that analysis points sit at least r_max from every boundary before computing K.

Antimeridian and projected-CRS distances. Computing inter-point distances in EPSG:4326 degrees near ±180° longitude wraps distances incorrectly, and Euclidean distance in degrees is not metric anyway. Project the pattern into an equal-area metric CRS such as EPSG:6933 (or a local UTM zone) before pair counting, so both the disc radius and the wrap arithmetic are in true metres.

Frequently Asked Questions

When should I use a guard region instead of an edge-corrected estimator?
Use a guard region when you generate the pattern yourself and can afford to synthesize a buffer wider than your largest search radius: analysing only interior points against the full buffered set removes the bias with no statistical correction and no assumptions about window shape. Use an edge-corrected estimator when the window is fixed and irregular, or when you cannot regenerate the pattern with a buffer.
Why does the bias get worse as the search radius grows?
The fraction of a point's search disc that falls outside the window grows with radius, so more of every boundary point's neighbourhood goes uncounted. Once the radius approaches the window's short dimension, a large share of every disc spills outside and the naive K collapses far below the theoretical value. That is why edge-correction tests should gate on the relative error at large radius.
Is toroidal wrapping ever wrong to use?
Yes. Toroidal wrapping assumes a stationary process on a rectangular window that tiles the plane. On an irregular window, a coastline-clipped area, or a pattern with an intensity gradient, it invents neighbour relationships across a boundary with no physical meaning and can bias the statistic upward. Restrict it to rectangular windows with verified stationarity.
Do edge effects matter if I only care about point count, not clustering?
First-order intensity — the expected number of points per unit area — is unbiased by the window boundary, so a simple count is fine. Edge effects only corrupt second-order statistics that count pairs or neighbour distances, such as Ripley's K and the nearest-neighbour distribution. If your realism gate scores clustering, you must correct for edges before comparing synthetic to source.