How to apply donut geomasking to sensitive points in Python

Problem statement

You have a point per subject — a case, an incident, a nest — and you need to publish points rather than counts. Donut geomasking displaces each one by a random distance between an inner and an outer radius, so no point stays where it was and none travels further than the analysis can tolerate.

Four details separate a working implementation from a broken one:

• The radius must be drawn with a square root, or points crowd towards the inner edge: with a uniform radius, 24.2% of points landed in the inner ninth of a 50–300 m ring's area, against 11.1% with the square root.
• The displacement must happen in a projected CRS, or the ring becomes an ellipse.
• The result must be constrained to plausible ground; an unconstrained 50–300 m donut put 77 of 4,064 Brighton & Hove address points into the sea.
• The output must be stored, not regenerated, because independent redraws average back to the truth.

This guide builds the function with all four handled, and measures what it did.

Quick answer

import numpy as np
import geopandas as gpd

def donut_mask(gdf, r_min, r_max, crs, seed):
    """Displace each point into the annulus [r_min, r_max], in metres."""
    rng = np.random.default_rng(seed)
    work = gdf.to_crs(crs)
    xy = np.c_[work.geometry.x, work.geometry.y]

    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)
    moved = xy + np.c_[r * np.cos(theta), r * np.sin(theta)]

    out = work.copy()
    out.geometry = gpd.points_from_xy(moved[:, 0], moved[:, 1], crs=crs)
    return out.to_crs(gdf.crs)

masked = donut_mask(cases, r_min=50, r_max=300, crs="EPSG:27700", seed=20260915)
d = cases.to_crs(27700).distance(masked.to_crs(27700))
print(f"displacement: median {d.median():.1f} m, min {d.min():.1f} m, max {d.max():.1f} m")

displacement: median 213.8 m, min 50.2 m, max 300.0 m

On the 4,064 Brighton & Hove address points this produced a median displacement of 213.8 m, a minimum of 50.2 m and a maximum of 300.0 m — the inner radius holding, which is the whole purpose of the donut.

Metric CRS, separate angle and radius, per-subject masks

1. Reproject to a metric CRS

Pick a projection appropriate to the extent — a national grid, a UTM zone, or an equidistant projection centred on the data. Masking in degrees stretches the ring by 1/cos(latitude): at 50.8°N a "300 m" east–west displacement is 475 m on the ground.

2. Draw the angle and radius separately

theta = rng.uniform(0, 2 * np.pi, n)
r = np.sqrt(rng.uniform(0, 1, n) * (r_max**2 - r_min**2) + r_min**2)

The square root converts a uniform draw on area into a uniform draw on radius. Without it, r = rng.uniform(r_min, r_max, n) concentrates points near r_min and the effective protection is much lower than the outer radius suggests.

3. Choose the radii from a candidate count

r_max should reach the number of candidate subjects you are promising. Compute the distance to the kth nearest dwelling and look at the distribution; on the Brighton points reaching k=5 needed a median of 14.0 m but a maximum of 421.0 m: even in this dense city a fixed 300 m leaves the most isolated address with fewer than five neighbours inside its outer radius. In a rural area the tail runs much further; use an adaptive one — step 7.

4. Mask per subject, not per record

Two events at the same address must move by the same vector. Group by subject, mask the group's location once, and apply it to every record.

subjects = cases.dissolve(by="subject_id").geometry.representative_point()
offsets = donut_offsets(len(subjects), 50, 300, rng)
lookup = dict(zip(subjects.index, offsets))

5. Reject points that land somewhere impossible

Water, airfields, the far side of an administrative boundary. Redraw only the offenders, and loop until none remain.

land = gpd.read_file("land.gpkg").geometry.union_all()
bad = ~masked.geometry.within(land)
rounds = 0
while bad.any() and rounds < 50:
    masked.loc[bad, "geometry"] = donut_mask(cases[bad], 50, 300, crs, seed + rounds).geometry
    bad = ~masked.geometry.within(land)
    rounds += 1

On the coastal Brighton set 77 points needed redrawing, one of them six times. The displacement distribution was barely affected: median 212.9 m, maximum 300.0 m.

6. Measure what the mask did

Report the displacement percentiles and at least one analysis-level effect. Geomasking changed the result of the analysis covers which statistics to check.

7. Consider the adaptive variant

Instead of one r_max for everybody, size each point's radius from its own kth-nearest-candidate distance. It moves the median point far less and still protects the isolated ones.

8. Store the output and the parameters

Write the masked file, and write a small JSON beside it with the method, radii, CRS, seed and date. Republishing means shipping the same file again.

Masking with a land constraint, end to end

The data used on this page

The points are every OpenStreetMap node in Brighton & Hove carrying an addr:housenumber — address points standing in for cases — and the candidate dwellings around them are every OpenStreetMap feature there with an address, by representative point. Land is the ONS December 2024 district boundaries clipped to the coastline (BFC) for Brighton & Hove and its neighbours, so a point that lands in the English Channel is rejected. Overpass via the maps.mail.ru mirror, read from its answer of 9 October 2026 (4,064 address points), which the geoprivacy guides share.

import geopandas as gpd
import io
import numpy as np
import os
import osmnx as ox
import requests
import zipfile
from scipy.spatial import cKDTree

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-brighton-addresses-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

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})
cases = features[features.index.get_level_values(0) == "node"][["geometry"]].reset_index(drop=True)
dwellings = features.to_crs(27700).representative_point()
dwellings_xy = np.c_[dwellings.x, dwellings.y]

w, s_, e, n_ = gpd.GeoSeries([area], crs=4326).to_crs(27700).buffer(2000).to_crs(4326).total_bounds
LAD = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Local_Authority_Districts_December_2024_Boundaries_UK_BFC/FeatureServer/0/query")
reply = requests.get(LAD, params={"geometry": f"{w},{s_},{e},{n_}", "geometryType": "esriGeometryEnvelope",
                                  "inSR": 4326, "outFields": "LAD24NM", "outSR": 27700, "f": "geojson"}, timeout=120)
land = gpd.read_file(reply.text)
print(f"{len(cases):,} address points, {len(dwellings_xy):,} addressed features; land from {sorted(land.LAD24NM)}")

4,064 address points, 61,142 addressed features; land from ['Adur', 'Brighton and Hove', 'Horsham', 'Lewes', 'Mid Sussex']

Example 1 — masking with a land constraint, end to end

import numpy as np, geopandas as gpd
from shapely.geometry import Point

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

def masked_on_land(gdf, land, r_min, r_max, crs="EPSG:27700", seed=0, max_rounds=50):
    rng = np.random.default_rng(seed)
    work, land = gdf.to_crs(crs), land.to_crs(crs).union_all()
    xy = np.c_[work.geometry.x, work.geometry.y]

    out = xy + donut_offsets(len(xy), r_min, r_max, rng)
    tries = np.ones(len(xy), int)
    bad = np.array([not land.contains(Point(p)) for p in out])
    for _ in range(max_rounds):
        if not bad.any():
            break
        idx = np.flatnonzero(bad)
        out[idx] = xy[idx] + donut_offsets(len(idx), r_min, r_max, rng)
        tries[idx] += 1
        bad[idx] = [not land.contains(Point(p)) for p in out[idx]]

    res = work.copy()
    res.geometry = gpd.points_from_xy(out[:, 0], out[:, 1], crs=crs)
    res["mask_attempts"] = tries
    return res.to_crs(gdf.crs)

masked = masked_on_land(cases, land, 50, 300, seed=20260915)
print(f"redrawn: {(masked.mask_attempts > 1).sum()}, max attempts {masked.mask_attempts.max()}")

redrawn: 77, max attempts 6

Example 2 — the adaptive radius

from scipy.spatial import cKDTree

def adaptive_donut(xy, candidates_xy, k=5, r_min=25, rng=None):
    rng = rng or np.random.default_rng(0)
    d, _ = cKDTree(candidates_xy).query(xy, k=k + 1)
    r_max = np.maximum(d[:, k], r_min * 2)             # never smaller than the inner ring
    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)]

xy = np.c_[cases.to_crs(27700).geometry.x, cases.to_crs(27700).geometry.y]
moved = adaptive_donut(xy, dwellings_xy, k=10, rng=np.random.default_rng(1))
d = np.hypot(*(moved - xy).T)
print(f"displacement: median {np.median(d):.1f} m, max {d.max():.1f} m")

displacement: median 40.5 m, max 261.1 m

Compare with 213.8 m median for the fixed 50–300 m donut: the adaptive version buys the same candidate count for a fifth of the median error, and moves its most isolated point 261 m.

Example 3 — recording the release

import json, hashlib, datetime, pathlib

masked.to_file("cases_masked.gpkg", layer="cases", driver="GPKG")

meta = {
    "method": "donut geomask with land-mask rejection",
    "r_min_m": 50, "r_max_m": 300,
    "crs": "EPSG:27700",
    "seed": 20260915,
    "records": int(len(masked)),
    "redrawn": int((masked.mask_attempts > 1).sum()),
    "displacement_median_m": round(float(cases.to_crs(27700).distance(masked.to_crs(27700)).median()), 1),
    "created": datetime.date.today().isoformat(),
    "sha256": hashlib.sha256(pathlib.Path("cases_masked.gpkg").read_bytes()).hexdigest(),
}
pathlib.Path("cases_masked.json").write_text(json.dumps(meta, indent=2))
print({k: meta[k] for k in ("records", "redrawn", "displacement_median_m")})

{'records': 4064, 'redrawn': 77, 'displacement_median_m': 212.9}

Example 4 — the radius and projection claims, measured

rng = np.random.default_rng(0)
n, r_min, r_max = 100_000, 50, 300
inner_ninth = np.sqrt(r_min**2 + (r_max**2 - r_min**2) / 9)        # radius enclosing the inner ninth of the ring's area
naive = rng.uniform(r_min, r_max, n)
right = np.sqrt(rng.uniform(0, 1, n) * (r_max**2 - r_min**2) + r_min**2)
print(f"inside the inner ninth of the ring's area: uniform radius {np.mean(naive < inner_ninth):.1%}, "
      f"square-root draw {np.mean(right < inner_ninth):.1%}")
print(f"for a full disc (r_min 0): uniform radius {np.mean(rng.uniform(0, r_max, n) < r_max / 3):.1%}")
lat = cases.geometry.y.mean()
print(f"at {lat:.1f}°N a '300 m' east-west step in degrees is {300 / np.cos(np.radians(lat)):.0f} m on the ground")
kth, _ = cKDTree(dwellings_xy).query(xy, k=6)
print(f"distance to the 5th nearest other address: median {np.median(kth[:, 5]):.1f} m, max {kth[:, 5].max():.1f} m")
for lo, hi in ((0, 100), (0, 300), (50, 300), (100, 500)):
    d = np.hypot(*donut_offsets(len(xy), lo, hi, np.random.default_rng(20260915)).T)
    print(f"{lo}-{hi} m: median {np.median(d):.1f} m, p5 {np.percentile(d, 5):.1f} m")

inside the inner ninth of the ring's area: uniform radius 24.2%, square-root draw 11.1%
for a full disc (r_min 0): uniform radius 33.3%
at 50.8°N a '300 m' east-west step in degrees is 475 m on the ground
distance to the 5th nearest other address: median 14.0 m, max 421.0 m
0-100 m: median 70.3 m, p5 21.2 m
0-300 m: median 210.8 m, p5 63.5 m
50-300 m: median 213.8 m, p5 80.2 m
100-500 m: median 358.4 m, p5 144.1 m

The checksum is what lets you prove next year that you republished the same file rather than a fresh draw.

The square root, rejection sampling and seeds

Why the square root is not optional

Points should be uniform over the area of the annulus. Area grows with r², so drawing r uniformly puts too many points near the inner edge. With r_min=0, r_max=300 and a linear draw, half the points land within 150 m; with the square-root draw, half land beyond 212 m. The measured median for the correct draw on these data was 211.3 m.

Why rejection sampling rather than snapping

Snapping an offending point to the nearest land moves it to the coastline, which is a visible, systematic artefact and tells an attacker the true point was inland. Redrawing keeps the distribution the mask promises. Three rounds sufficed for a coastal city; if a dataset needs dozens, the radius is too large for the geography.

Why the seed belongs in the metadata

A seed is not a secret. Publishing it lets a user reproduce your file exactly and lets you prove that the version you shipped in two consecutive years is the same draw. What must not happen is a new draw each year: sixteen independent 50–300 m masks averaged to within 50 m of the truth for 59.0% of points.

Why the mask cannot fix a dataset that is too small

Masking redistributes uncertainty; it does not create candidates. Twenty cases in a village of thirty houses stay identifiable at any radius that keeps them in the village. In that situation aggregate, or publish the attribute at a coarser geography.

Redraw caps, derived data and coordinate grids

• Cap the redraw loop. An impossible constraint should raise, not spin.
• Mask before you derive anything. Centroids, densities and counts must come from the masked points.
• Keep the CRS in the metadata. The radii are meaningless without it.
• Do not round the masked coordinates to a grid. It reintroduces a visible lattice.
• Multi-point subjects move together. One offset per subject, applied to every record.
• Check the bounding box afterwards. A mask enlarges the extent by r_max on every side.
• Points near a boundary can cross it. If a region column is also published, constrain the draw to the region.
• A masked release is still personal data. Treat the file accordingly.

Internal links

• Geomasking methods explained: donut, random and adaptive — the family and the parameter choices
• Geomasked points land in the sea or outside the study area — when the rejection loop is not enough
• Spatial k-anonymity explained — choosing the outer radius
• Geomasking changed the result of the analysis — reporting the cost
• How to measure re-identification risk in a point dataset — the audit before and after
• How to choose a projected CRS for an area — working in metres
• How to reproject spatial data in GeoPandas — the round trip
• How to run a privacy check before publishing a spatial dataset — the gate that checks the mask ran

FAQ

What inner and outer radius should I use?

The outer radius from the candidate count you are promising, the inner from the size of a parcel — 25 m in dense terraces, 50–100 m in suburbs. Publish both.

Why do my masked points bunch near the inner edge?

Because the radius was drawn uniformly. Use r = sqrt(U*(r_max² − r_min²) + r_min²) so points are uniform over the area of the ring.

Can I mask in EPSG:4326?

No. Adding degrees produces an ellipse that is wrong by 1/cos(latitude) east–west. Reproject to a metric CRS, mask, reproject back.

What do I do about points that land in water?

Redraw them. On the Brighton address points this affected 77 of 4,064, with barely any change to the displacement distribution.

Should I re-run the mask for the next release?

No. Publish the same masked file. Independent redraws let an attacker average towards the truth — 59% of points within 50 m after sixteen releases in the test above.

Is donut masking better than aggregation?

Only if the analysis genuinely needs points. Aggregation removes the individual record, which no masking method does.