Geomasking changed the result of the analysis

Problem statement

The masked dataset went out, somebody ran the analysis they had run on the unmasked version, and the answer was different. The hotspots moved. The nearest-neighbour statistic said the points were no longer clustered. The counts per ward did not match the counts you published.

This is not a bug. Masking adds a known amount of random error to every coordinate, and every spatial statistic is a function of those coordinates. What is a bug is shipping the masked data without saying how much it moved, so that users treat masked results as unmasked ones.

On 4,059 real address points in Brighton & Hove, a 50–300 m donut mask raised the mean nearest-neighbour distance from 13.7 m to 36.6 m, left only 68.8% of the top-5% kernel density cells in the top 5%, and moved 45.6% of points into a different 500 m cell.

Quick answer

Measure the distortion for the statistics your users will compute, and publish the table:

import numpy as np
from scipy.spatial import cKDTree
from scipy.stats import spearmanr

def distortion(true_xy, masked_xy, grid, bw=0.05):
    nn_t = cKDTree(true_xy).query(true_xy, k=2)[0][:, 1]
    nn_m = cKDTree(masked_xy).query(masked_xy, k=2)[0][:, 1]
    kde_t, kde_m = kde(true_xy, grid, bw), kde(masked_xy, grid, bw)
    top = kde_t >= np.percentile(kde_t, 95)
    return {
        "mean_nn_ratio": nn_m.mean() / nn_t.mean(),
        "kde_spearman": spearmanr(kde_t, kde_m).statistic,
        "top5_retained": float(((kde_m >= np.percentile(kde_m, 95)) & top).sum() / top.sum()),
    }

Publish it next to the data. A user who knows the nearest-neighbour statistic is inflated 2.7× can correct for it; a user who does not will report the inflation as a finding.

Measure damage across radii, publish true aggregates

1. Measure the damage across the radius range

The full curve on the Brighton address points:

Even the smallest mask nearly doubles the nearest-neighbour statistic. The coarse-bandwidth density survives far better than the fine one, which is the clue to the mitigation.

2. Check the aggregate counts separately

Masking moves points across cell boundaries, so any count computed from the masked points differs from the true one:

3. Publish the aggregates yourself

The fix for the count mismatch is not a better mask. It is to compute the aggregates from the true points, publish those, and mask only the point layer — so users take counts from the table and pattern from the points.

4. Constrain the mask to the reporting unit

If the release includes a region column, masking within the region keeps every published total exact. The cost is that the region is disclosed, which it already was.

5. Match the analysis scale to the mask

A mask of radius R destroys structure below roughly R and preserves structure well above it. Tell users the resolution the data supports: "do not interpret features smaller than 500 m". The coarse-bandwidth density above kept ρ=0.9945 under the same mask that took the fine one to 0.9775 — and the difference is sharper in the hotspot cells, where the fine surface kept only 68.8% of its top 5%.

6. Use an adaptive radius where accuracy matters most

An adaptive mask spends displacement where candidates are scarce, so dense areas — usually the interesting ones — barely move. On these points the median displacement at k=10 was 40.5 m against 213.8 m for the fixed 50–300 m donut (the donut guide has the code).

7. Correct the statistics where a correction exists

For some statistics the bias is analytic. Masking adds independent noise of known variance, so second-moment measures can be deconvolved; a K-function or a variogram can have the noise variance subtracted. Publish the noise variance so users can do it.

A distortion sweep on Brighton addresses

The data used on this page

Every OpenStreetMap node in Brighton & Hove carrying an addr:housenumber — 4,059 address points standing in for cases — in British National Grid metres, and the ONS December 2024 wards of Brighton and Hove for the per-ward examples. Each mask is one draw with a fixed seed. Overpass via the maps.mail.ru mirror.

import geopandas as gpd
import numpy as np
import osmnx as ox
import pandas as pd
import requests
from shapely import contains, points as spoints

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.overpass_rate_limit = False
ox.settings.http_user_agent = "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"
area = ox.geocode_to_gdf("Brighton and Hove, England").geometry.iloc[0]
features = ox.features_from_polygon(area, tags={"addr:housenumber": True})
true_points = features[features.index.get_level_values(0) == "node"][["geometry"]].reset_index(drop=True).to_crs(27700)
xy = np.c_[true_points.geometry.x, true_points.geometry.y]
rng = np.random.default_rng(20260915)


def donut(xy, r_min, r_max):
    theta = rng.uniform(0, 2 * np.pi, len(xy))
    r = np.sqrt(rng.uniform(0, 1, len(xy)) * (r_max**2 - r_min**2) + r_min**2)
    return xy + np.c_[r * np.cos(theta), r * np.sin(theta)]


WARDS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
         "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query")
wards = gpd.read_file(requests.get(WARDS, params={"where": "LAD24NM='Brighton and Hove'", "outFields": "WD24NM",
                                                  "outSR": 27700, "f": "geojson"}, timeout=120).text)
wards = wards.rename(columns={"WD24NM": "ward"})[["ward", "geometry"]]


def mask_within_regions(points, regions, r_min, r_max, max_rounds=100):
    """Donut-mask each point, redrawing until it lands in the region it started in."""
    home = gpd.sjoin(points, regions, predicate="within")["ward"]
    pts = points.loc[home.index]
    start = np.c_[pts.geometry.x, pts.geometry.y]
    shapes = regions.set_index("ward").geometry.reindex(home.values).values
    out = donut(start, r_min, r_max)
    bad = ~contains(shapes, spoints(out))
    for _ in range(max_rounds):
        if not bad.any():
            break
        out[bad] = donut(start[bad], r_min, r_max)
        bad = ~contains(shapes, spoints(out))
    return gpd.GeoDataFrame(geometry=gpd.points_from_xy(out[:, 0], out[:, 1]), crs=points.crs, index=pts.index)


print(f"{len(true_points):,} address points, {len(wards)} wards")

4,059 address points, 23 wards

Example 1 — the full distortion sweep

import numpy as np
from scipy.spatial import cKDTree
from scipy.stats import gaussian_kde, spearmanr

gx, gy = np.meshgrid(np.linspace(xy[:, 0].min(), xy[:, 0].max(), 120),
                     np.linspace(xy[:, 1].min(), xy[:, 1].max(), 120))
grid = np.vstack([gx.ravel(), gy.ravel()])

def kde(points, bw):
    return gaussian_kde(points.T, bw_method=bw)(grid)

truth = {bw: kde(xy, bw) for bw in (0.05, 0.15)}
base_nn = cKDTree(xy).query(xy, k=2)[0][:, 1]

top_true = truth[0.05] >= np.percentile(truth[0.05], 95)
cell_true = pd.Series(list(zip((xy[:, 0] // 500).astype(int), (xy[:, 1] // 500).astype(int))))
results = {}
print(f"none            {base_nn.mean():6.1f}")
for label, (r0, r1) in [("donut 0-100", (0, 100)), ("donut 50-300", (50, 300)),
                        ("donut 100-500", (100, 500)), ("donut 200-1000", (200, 1000))]:
    m = donut(xy, r0, r1)
    nn = cKDTree(m).query(m, k=2)[0][:, 1]
    row = [f"{label:14}", f"{nn.mean():6.1f}", f"{nn.mean() / base_nn.mean():6.3f}"]
    rho = [spearmanr(truth[bw], kde(m, bw)).statistic for bw in (0.05, 0.15)]
    row += [f"{r:8.4f}" for r in rho]
    k_m = kde(m, 0.05)
    kept = ((k_m >= np.percentile(k_m, 95)) & top_true).sum() / top_true.sum()
    cell_m = pd.Series(list(zip((m[:, 0] // 500).astype(int), (m[:, 1] // 500).astype(int))))
    ct, cm = cell_true.value_counts(), cell_m.value_counts()
    idx = ct.index.union(cm.index)
    mae = (ct.reindex(idx, fill_value=0) - cm.reindex(idx, fill_value=0)).abs().mean()
    moved = (cell_true != cell_m).mean()
    row += [f"top5 kept {kept:6.1%}", f"cells {len(idx)}", f"MAE {mae:5.1f}", f"moved {moved:6.1%}"]
    results[label] = dict(ratio=nn.mean() / base_nn.mean(), rho_fine=rho[0], rho_coarse=rho[1], kept=kept, moved=moved)
    print("  ".join(row))

none              13.7
donut 0-100       25.5   1.868    0.9970    0.9995  top5 kept  91.5%  cells 160  MAE   2.8  moved  16.3%
donut 50-300      36.6   2.678    0.9775    0.9945  top5 kept  68.8%  cells 186  MAE   7.1  moved  45.6%
donut 100-500     42.8   3.134    0.9501    0.9804  top5 kept  54.3%  cells 213  MAE  10.6  moved  70.5%
donut 200-1000    53.1   3.883    0.8536    0.9020  top5 kept  44.2%  cells 275  MAE  13.5  moved  94.5%

Example 2 — counts from the truth, points from the mask

import geopandas as gpd

# Aggregate the true points and publish that table
true_counts = gpd.sjoin(true_points, wards, predicate="within") \
                 .groupby("ward").size().rename("cases")

# Mask within each ward so the masked points also fall in the right one
masked = mask_within_regions(true_points, wards, r_min=50, r_max=300)
masked_counts = gpd.sjoin(masked, wards, predicate="within").groupby("ward").size()

assert (true_counts == masked_counts.reindex(true_counts.index, fill_value=0)).all()
print(f"{len(true_counts)} wards, {true_counts.sum():,} points: published and masked counts agree")

23 wards, 4,059 points: published and masked counts agree

The assertion is the point: if masking is constrained to the reporting units, the published table and the published points agree exactly, and users stop filing bug reports.

Example 3 — the statement to ship with the data

import json

statement = {
    "method": "donut geomask, land-constrained, per-ward",
    "r_min_m": 50, "r_max_m": 300,
    "effect": {
        "mean_nearest_neighbour_ratio": round(results["donut 50-300"]["ratio"], 3),
        "kde_spearman_fine_bandwidth": round(results["donut 50-300"]["rho_fine"], 4),
        "kde_spearman_coarse_bandwidth": round(results["donut 50-300"]["rho_coarse"], 4),
        "top5_density_cells_retained": round(results["donut 50-300"]["kept"], 3),
        "points_in_a_different_500m_cell": round(results["donut 50-300"]["moved"], 3),
    },
    "guidance": "Do not interpret features smaller than about 500 m. "
                "Ward totals are published separately and are exact.",
}
print(json.dumps(statement["effect"], indent=2))

{
  "mean_nearest_neighbour_ratio": 2.678,
  "kde_spearman_fine_bandwidth": 0.9775,
  "kde_spearman_coarse_bandwidth": 0.9945,
  "top5_density_cells_retained": 0.688,
  "points_in_a_different_500m_cell": 0.456
}

Why nearest-neighbour statistics break first

Why nearest-neighbour statistics break first

Masking adds error of the order of the radius, so any statistic computed at a scale smaller than the radius is dominated by it. The mean nearest-neighbour distance on these addresses was 13.7 m — some fifteen times smaller than the mask's median displacement — so the measured value afterwards is essentially the mask, not the data. Anything based on short distances (clustering indices, duplicate detection, building-level joins) is in the same position.

Why coarse density survives and fine density does not

A kernel density estimate at bandwidth h is a smoothed version of the point pattern. Masking with radius R convolves the pattern with a second kernel of roughly that width. When h ≫ R, the extra convolution barely changes the result: ρ = 0.9945 at the coarse bandwidth. When h ≈ R, the two are comparable and the surface changes more: ρ = 0.9775 overall, and a third of the hotspot cells lost. The rule of thumb is that the smallest interpretable feature is about the mask radius.

Why the counts must come from the unmasked data

There is no masking radius small enough to keep points in their original cells and large enough to protect them — 9.1% of points crossed a 500 m boundary under a 100 m mask. Either constrain the mask to the units, or publish the counts separately. Doing neither guarantees that the two published products contradict each other.

Why the statement matters as much as the method

A masked dataset with no accuracy statement will be used as if it were exact, and the first result anyone publishes from it will be wrong in a way that is attributable to you. The table above takes an afternoon to compute and removes that entire class of problem.

Hotspots, distance analyses and interpolation

• Hotspot tests inflate their false-positive rate. Masking flattens peaks; significance thresholds calibrated on unmasked data no longer hold.
• Distance-to-feature analyses are badly hit. "Within 100 m of a road" is meaningless under a 300 m mask.
• Point-in-polygon joins move. Any attribute enriched after masking is enriched from the wrong polygon.
• Interpolation smooths twice. A surface built from masked points is smoother than its stated bandwidth.
• Pair counts are inflated at short range. Ripley's K is biased downward near zero.
• Small samples are affected more. With few points, one displacement changes the statistic.
• The mask variance is publishable. For a donut it is a closed form, and users can deconvolve with it.
• Aggregations coarser than the radius are safe. That is the one thing you can promise without caveat.

Internal links

• Geomasking methods explained: donut, random and adaptive — choosing the radius with this cost in mind
• How to apply donut geomasking to sensitive points in Python — the implementation, including per-region masking
• Spatial k-anonymity explained — the adaptive radius that costs far less accuracy
• Kernel density explained — bandwidth and what it can resolve
• KDE bandwidth is wrong — choosing h for masked data
• The modifiable areal unit problem explained — the other reason counts move
• How to find hotspots with Getis-Ord in Python — the test whose false-positive rate changes
• Geoprivacy explained: why coordinates are personal data — when to aggregate instead of masking

FAQ

How much does geomasking change my results?

It depends on the radius and the statistic. A 50–300 m mask on real address points raised the mean nearest-neighbour distance 2.7×, kept 68.8% of the top-5% density cells and moved 45.6% of points into a different 500 m cell.

Why do my masked counts not match the published totals?

Because masking moves points across boundaries. Publish counts computed from the true points, and constrain the mask to the reporting units if the points must also aggregate correctly.

Which analyses survive masking?

Anything computed at a scale well above the radius. Coarse density kept a rank correlation of 0.9945 under a mask that took the fine-bandwidth version to 0.9775 and its hotspots to 68.8%.

Can I correct for the bias?

Partly. The mask adds independent noise of known variance, so second-moment statistics such as variograms and K-functions can be corrected if you publish the variance.

Is there a mask that does not distort anything?

No. Every displacement changes every statistic. What you can choose is where the error goes — an adaptive radius puts it in sparse areas and leaves dense ones nearly intact.

What should I publish alongside masked points?

The method and parameters, the measured distortion for the common statistics, the exact aggregate counts, and a plain statement of the smallest feature users should interpret.