Correcting Spatial Sampling Bias with GeoPandas

TL;DR: Call scipy.spatial.Voronoi on your projected point coordinates, clip each cell to the study boundary with geopandas.GeoSeries.intersection, compute weight = 1 / cell_area, normalise to mean = 1.0, then clip to [0.1, 10.0]. The resulting spatial_weight column passes directly into statsmodels WLS, XGBoost’s sample_weight, or a bias-corrected subsample for variogram fitting. The full function is in Step 3 below.

Why This Matters #

Non-uniform data collection creates systematic overrepresentation of accessible areas — road corridors, urban cores, well-funded monitoring stations — while remote or ecologically sensitive zones are undersampled. Feeding that raw point cloud into ordinary or universal kriging, spatial regression models, or ecological niche algorithms produces biased parameter estimates and inflated prediction confidence in the oversampled zones. Sampling bias mitigation addresses this before any model is fit, ensuring that each observation’s statistical contribution reflects its spatial uniqueness rather than its proximity to the nearest field crew. This page shows the production-grade implementation within the broader core concepts of spatial statistics and geostatistics framework.

Environment and Version Pinning #

# pip install geopandas==0.14.4 scipy==1.13.0 shapely==2.0.4 numpy==1.26.4
import numpy as np                          # 1.26.4
import geopandas as gpd                     # 0.14.4
from scipy.spatial import Voronoi           # 1.13.0
from shapely.geometry import Polygon        # 2.0.4
from shapely.ops import unary_union

CRS note: every snippet below assumes a projected CRS in metres (e.g., UTM or an equal-area projection). Geographic coordinates (EPSG:4326) will silently produce degree-squared areas — always reproject first.

How Voronoi Inverse-Probability Weighting Works #

Before writing code it helps to see the geometric logic. Each sample point “owns” the region of space closer to it than to any other point. That region’s area is the inverse of local sampling intensity: a tiny cell means the point is in a densely sampled neighbourhood; a large cell means it stands alone. Weighting by the reciprocal of area therefore gives isolated points more influence and dense clusters less.

The mathematical formula for the weight assigned to point $i$ is:

$$w_i = \frac{1}{A_i}, \quad \tilde{w}_i = \frac{w_i}{\bar{w}}, \quad \hat{w}_i = \text{clip}\!\left(\tilde{w}_i,\, w_{\min},\, w_{\max}\right)$$

where $A_i$ is the clipped Voronoi cell area for point $i$, $\bar{w}$ is the mean raw weight across all points, and the final clipped weight $\hat{w}_i$ is bounded to $[w_{\min}, w_{\max}]$ for numerical stability.

Step-by-Step Implementation #

Step 1 — Load and Reproject Data #

import geopandas as gpd

# Load your point observations and study boundary
points = gpd.read_file("observations.gpkg")
boundary = gpd.read_file("study_area.gpkg")

# Reproject to a projected CRS; UTM zone 32N as an example
# Replace EPSG:32632 with the UTM zone covering your study area
target_crs = "EPSG:32632"
points = points.to_crs(target_crs)
boundary = boundary.to_crs(target_crs)

# Quick sanity check
assert not points.crs.is_geographic, "CRS must be projected (metres)."
print(f"Points: {len(points)} observations, CRS: {points.crs.to_epsg()}")

Choosing the right projection matters: equal-area projections (e.g., EPSG:6933 NSIDC EASE-Grid 2.0) minimise area distortion across large extents, while UTM zones are more accurate within a single 6° strip.

Step 2 — Build the Voronoi Tessellation #

import numpy as np
from scipy.spatial import Voronoi
from shapely.geometry import Polygon
import geopandas as gpd

def _build_voronoi_gdf(points_gdf, boundary_geom):
    """
    Returns a GeoDataFrame of clipped Voronoi polygons,
    index-aligned with points_gdf.
    """
    coords = np.column_stack([
        points_gdf.geometry.x.values,
        points_gdf.geometry.y.values
    ])
    n = len(coords)
    if n < 3:
        raise ValueError("At least 3 non-collinear points required.")

    vor = Voronoi(coords)  # O(n log n)

    # Bounding box to capture infinite edge cells
    minx, miny, maxx, maxy = boundary_geom.bounds
    pad = max(maxx - minx, maxy - miny)
    bbox = Polygon([
        (minx - pad, miny - pad), (maxx + pad, miny - pad),
        (maxx + pad, maxy + pad), (minx - pad, maxy + pad)
    ])

    polys = {}
    for pt_idx, reg_idx in enumerate(vor.point_region):
        region = vor.regions[reg_idx]
        if not region or -1 in region:
            # Edge cell: use bbox as conservative bounds
            polys[pt_idx] = bbox
        else:
            polys[pt_idx] = Polygon([vor.vertices[v] for v in region])

    geom_series = [polys[i] for i in range(n)]
    vor_gdf = gpd.GeoDataFrame(geometry=geom_series,
                                index=range(n),
                                crs=points_gdf.crs)
    # Clip every cell to the study boundary
    vor_gdf["geometry"] = vor_gdf.geometry.intersection(boundary_geom)
    return vor_gdf

vor.point_region is a one-to-one mapping from input point index to region index — it ensures the clipped polygon list stays aligned with the original point order.

Step 3 — Compute and Normalise Weights #

def compute_voronoi_weights(points_gdf, boundary_gdf,
                             min_weight=0.1, max_weight=10.0):
    """
    Adds a 'spatial_weight' column to points_gdf.

    Parameters
    ----------
    points_gdf : gpd.GeoDataFrame
        Projected point observations.
    boundary_gdf : gpd.GeoDataFrame
        Study area polygon(s); dissolved to a single geometry internally.
    min_weight, max_weight : float
        Stability bounds applied after normalisation.

    Returns
    -------
    gpd.GeoDataFrame with original columns plus 'spatial_weight'.
    """
    if points_gdf.crs.is_geographic:
        raise ValueError("Reproject to a metric CRS before calling this function.")

    # Dissolve to one polygon; buffer(0) heals any topology issues
    boundary_geom = boundary_gdf.geometry.union_all()
    if not boundary_geom.is_valid:
        boundary_geom = boundary_geom.buffer(0)

    vor_gdf = _build_voronoi_gdf(points_gdf, boundary_geom)

    # Area in CRS units² (metres² for UTM)
    vor_gdf["area_m2"] = vor_gdf.geometry.area
    vor_gdf["area_m2"] = vor_gdf["area_m2"].replace(0.0, np.nan)  # guard against degenerate cells

    # Inverse-area weight, normalised to mean = 1.0
    vor_gdf["weight"] = 1.0 / vor_gdf["area_m2"]
    mean_w = vor_gdf["weight"].mean(skipna=True)
    vor_gdf["weight"] = vor_gdf["weight"] / mean_w

    # Stability clipping
    vor_gdf["weight"] = vor_gdf["weight"].clip(lower=min_weight, upper=max_weight)

    # Map weights back to original GeoDataFrame
    result = points_gdf.copy()
    result["spatial_weight"] = vor_gdf["weight"].reindex(result.index,
                                                          fill_value=1.0)
    return result

Step 4 — Validate the Weight Distribution #

weighted_pts = compute_voronoi_weights(points, boundary)

# Summary statistics
print(weighted_pts["spatial_weight"].describe())

# Effective sample size — how much unique information the weighted sample carries
n = len(weighted_pts)
w = weighted_pts["spatial_weight"].values
ess = n / ((w ** 2).sum() / n)
print(f"ESS: {ess:.1f}  /  {n} raw observations  ({100*ess/n:.0f} %)")

Interpreting the Output #

A weight above 1.0 means the point is more spatially isolated than average and should count more in downstream fitting. A weight below 1.0 means it sits in a dense neighbourhood. The normalisation guarantees that the sum of weights equals $n$, preserving the scale of regression coefficients.

Critical Best Practices #

Always Project Before Area Calculations #

Geographic coordinates produce degree-squared areas — essentially meaningless for weighting because the physical size of a degree varies with latitude. assert not points_gdf.crs.is_geographic at function entry catches this early.

Guard Against Duplicate Coordinates #

Two points at identical locations produce a Voronoi cell of area zero, making the inverse-area weight infinite. Deduplicate or add a small spatial jitter (points.geometry = points.geometry.apply(lambda g: g.buffer(0.01).centroid)) before running tessellation.

Cap Extreme Weights #

Without clip(lower=0.1, upper=10.0), a single isolated point near the study boundary can receive a weight several hundred times the mean, dominating variance estimation in kriging or inflating regression coefficients. The cap values are problem-specific; for ecological data with strong spatial gradients, max_weight=5.0 is common.

Align Indices Before Merging #

vor.point_region guarantees positional alignment, not label alignment. Use .reindex(result.index) — not .values — to safely map weights back when the input GeoDataFrame has a non-default integer index.

Validate Against Cross-Validation Strategies #

After weighting, run a spatial $k$-fold split and compare weighted vs unweighted RMSE on held-out blocks. If the weighted model degrades on the test blocks, the sampling density pattern is actually informative (e.g., higher measurement effort in high-variability zones) and inverse-probability weighting would remove a legitimate signal.

Applying Weights to Downstream Models #

Once spatial_weight is attached to your GeoDataFrame, integration into common model fitting APIs is straightforward:

import statsmodels.formula.api as smf

# Weighted least squares spatial regression
wls_model = smf.wls(
    "log_pollution ~ elevation + distance_to_road",
    data=weighted_pts,
    weights=weighted_pts["spatial_weight"]
).fit()
print(wls_model.summary())

from xgboost import XGBRegressor

X = weighted_pts[["elevation", "distance_to_road"]].values
y = weighted_pts["log_pollution"].values
w = weighted_pts["spatial_weight"].values

xgb = XGBRegressor(n_estimators=400, learning_rate=0.05)
xgb.fit(X, y, sample_weight=w)

For stationarity testing and variogram fitting, pass a bias-corrected subsample: draw points with np.random.choice(n, size=int(ess), replace=False, p=w/w.sum()) to produce a spatially representative input without artificially inflating the point count.

Troubleshooting #

Next Steps #

Return to Sampling Bias Mitigation for the full range of correction strategies including kernel density and stratified resampling. For a worked example combining bias-corrected samples with kriging uncertainty surfaces, see Step-by-Step Ordinary Kriging with PyKrige.

Related #

• Sampling Bias Mitigation — parent page covering the full spectrum of bias correction strategies
• Ripley’s K-Function Implementation Guide — quantifying clustering intensity before deciding on a correction strategy
• Spatial K-Fold Cross-Validation Setup — validating that bias-corrected models generalise across space

← Back to Sampling Bias Mitigation