⚙️ How-To Privacy & sensitive locations Geometry & topology Geoprivacy & Anonymisation · step 8 of 20 Updated

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 a third of the points land in the inner ninth of the ring.
  • 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 63 of 3,109 Brighton points into the English Channel.
  • 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)

On the 3,109 Brighton address points this produced a median displacement of 215.0 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.

Flow from points through reprojection, annulus draw, land-mask rejection and storage of the masked file.
Reject-and-redraw sits between the draw and the write; nothing after it is random.

Step-by-step solution

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 21.9 m and a maximum of 533.6 m, so a fixed 300 m under-protects the isolated tail. If the tail matters, use an adaptive radius — 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 this needed 3 rounds, with a maximum of 4 attempts for one point and 63 points redrawn at all. The displacement distribution was unaffected: median 211.4 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.

Bars of median displacement and fifth-percentile displacement for four donut radius settings.
The fifth percentile is the protection the weakest points get; the median is what the analysis pays.

Code examples

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: 63, max attempts 4

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)]

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 23.6 m, max 580.1 m

Compare with 215.0 m median for the fixed 50–300 m donut: the adaptive version buys the same protection for a tenth of the median error.

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": 211.4,
    "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))

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

Explanation

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.

Triage of four donut masking mistakes — a linear radius draw, masking in degrees, one offset per record, and regenerating the mask — with the fix for each.
Every one of these produces a file whose displacement summary looks right.

Edge cases or notes

  • 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.

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 a coastal dataset this affected 63 of 3,109 points and resolved in three rounds, with no 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.