Geomasked points land in the sea or outside the study area #

Problem statement #

The mask ran, the displacement statistics look right, and the map shows a scatter of cases floating in the English Channel. On a coastal address dataset, a 50–300 m donut mask put 73 of 3,109 points (2.35%) into the water on one draw, and a 200–1000 m mask put 544 (17.5%) there.

This is not only ugly. A point in the sea is a point an attacker can rule out, and knowing the true location is on the landward side of the released one shrinks the candidate set — so the visible failure is also a real loss of protection. The same applies to points that cross an administrative boundary whose counts you are also publishing: the totals stop matching, and the mismatch localises the error.

Quick answer #

Reject and redraw against a land mask rather than snapping to the nearest land:

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

land_geom = land.to_crs(27700).union_all()

out = xy + donut_offsets(len(xy), 50, 300, rng)
bad = np.array([not land_geom.contains(Point(p)) for p in out])
rounds = 0
while bad.any() and rounds < 50:
    idx = np.flatnonzero(bad)
    out[idx] = xy[idx] + donut_offsets(len(idx), 50, 300, rng)
    bad[idx] = [not land_geom.contains(Point(p)) for p in out[idx]]
    rounds += 1

On the Brighton set this converged in 3 rounds, redrew 63 of 3,109 points, and needed at most 4 attempts for a single point. The displacement distribution was unchanged: median 211.4 m, maximum 300.0 m.

Step-by-step solution #

1. Confirm how many points are affected, and where #

in_water = masked.geometry.within(sea_geom)
print(f"{in_water.sum()} of {len(masked)} ({in_water.mean():.2%}) in water")

If the share is small and clustered along the coast, it is the radius meeting the geography. If points are in water everywhere, including inland, the mask is being applied in the wrong CRS.

2. Check the CRS before anything else #

Masking in degrees adds a number meant as metres to a coordinate in degrees. A "300 m" displacement then moves the point 300 degrees, which lands somewhere that is not on the planet, or — with a scaled radius — 1/cos(latitude) too far east–west. Reproject to a metric CRS, mask, reproject back.

3. Measure the coastal band #

The share of points at risk is the share close enough to water for the outer radius to reach it. On the Brighton addresses, 739 of 3,109 (23.8%) lay within 300 m of the sea, and roughly a tenth of those ended up in the water on a 50–300 m draw:

4. Reject and redraw, do not snap #

Snapping an offending point to the nearest land puts it on the coastline. Every masked point from the seafront then sits on one line, which is a visible artefact and a strong hint about the true position. Redrawing preserves the distribution the mask promises.

5. Cap the loop and fail loudly #

An impossible constraint — a mask radius that cannot fit in the geometry — must raise, not spin. A cap of 50 rounds is generous; if it is reached, the radius is wrong for the geography.

6. Use a land mask at the right resolution #

Natural Earth 10m coastline is a global dataset generalised for world maps and is tens to hundreds of metres out at city scale. Use a national coastline, an OpenStreetMap coastline extract, or a land-use polygon set for anything where the shoreline matters.

7. Constrain to the region as well, when regions are published #

If the release also carries a ward or district column, the mask must stay inside the region or the counts and the points disagree. Intersect the constraint geometry with the region for each point.

8. Consider an adaptive radius instead #

The points at risk are the ones near water, and they are often in dense areas where a small radius already reaches enough candidates. An adaptive mask sized from local density moves those points far less and avoids most of the problem — median displacement 23.6 m at k=10, against 215.0 m for the fixed donut.

Code examples #

Example 1 — build a land mask from an OpenStreetMap coastline #

import json, geopandas as gpd
from shapely.geometry import LineString, Polygon
from shapely.ops import linemerge, unary_union

ways = json.loads(open("coastline.json").read())["elements"]
lines = [LineString([(p["lon"], p["lat"]) for p in w["geometry"]])
         for w in ways if w.get("geometry")]
merged = linemerge(unary_union(lines))
coast = max(merged.geoms, key=lambda g: g.length) if merged.geom_type != "LineString" else merged

# close the coastline against the far edge of the bounding box to make a sea polygon
xs = [c[0] for c in coast.coords]
sea = Polygon(list(coast.coords) + [(xs[-1], BOX_SOUTH), (xs[0], BOX_SOUTH)])
sea = gpd.GeoDataFrame(geometry=[sea], crs=4326).to_crs(27700)
print(f"sea polygon: {sea.area.iloc[0] / 1e6:.1f} km², valid={sea.is_valid.iloc[0]}")

sea polygon: 167.0 km², valid=True

OpenStreetMap coastline ways run with land on the left, so closing the merged line against the seaward edge of the box gives the water polygon directly.

Example 2 — reject-and-redraw with an attempt counter #

import numpy as np
from shapely.geometry import Point
from shapely.prepared import prep

def mask_within(xy, offsets_fn, allowed, max_rounds=50):
    allowed = prep(allowed)                       # prepared geometry: ~10x faster contains
    out = xy + offsets_fn(len(xy))
    tries = np.ones(len(xy), int)
    bad = np.array([not allowed.contains(Point(p)) for p in out])
    for _ in range(max_rounds):
        if not bad.any():
            return out, tries
        idx = np.flatnonzero(bad)
        out[idx] = xy[idx] + offsets_fn(len(idx))
        tries[idx] += 1
        bad[idx] = [not allowed.contains(Point(p)) for p in out[idx]]
    raise RuntimeError(f"{bad.sum()} point(s) could not be placed in {max_rounds} rounds")

moved, tries = mask_within(xy, lambda n: donut_offsets(n, 50, 300, rng), land_geom)
print(f"redrawn: {(tries > 1).sum()}, max attempts {tries.max()}")

redrawn: 63, max attempts 4

shapely.prepared.prep matters here: the loop tests thousands of points against one large polygon, which is exactly the case prepared geometries exist for.

Example 3 — constrain to the region as well #

import geopandas as gpd

regions = gpd.read_file("wards.gpkg").to_crs(27700)
joined = gpd.sjoin(points, regions[["ward", "geometry"]], predicate="within")

moved = points.geometry.values.copy()
for ward, group in joined.groupby("ward"):
    allowed = regions.loc[regions.ward == ward, "geometry"].iloc[0].intersection(land_geom)
    xy = np.c_[group.geometry.x, group.geometry.y]
    out, _ = mask_within(xy, lambda n: donut_offsets(n, 50, 300, rng), allowed)
    moved[group.index] = gpd.points_from_xy(out[:, 0], out[:, 1])

Masking per region keeps every published aggregate consistent with the published points, at the cost of telling an attacker which region each point is in — which the release was publishing anyway.

Explanation #

Why snapping is worse than redrawing #

Rejection sampling produces exactly the intended distribution restricted to the allowed area. Snapping produces a spike on the boundary: every rejected point ends up on the coastline, a shape nobody would mistake for data, and the distance from the coast to the true point is bounded by the radius. The visible artefact and the information leak are the same defect.

Why the failures cluster and that is fine #

Only points within r_max of water can land in it, so on this dataset 23.8% of points were exposed at a 300 m radius and about a tenth of those failed. The redraw cost is therefore small and bounded; if it is not, the radius is too large for the coastline and the mask needs rethinking rather than more rounds.

Why the coastline's resolution matters #

A generalised global coastline can sit a hundred metres from the real one, so points will be rejected on land and accepted in water. If the shoreline is where the data lives — coastal towns, ports, beaches — use a coastline whose accuracy is better than the mask radius.

Why the same problem appears with lakes, airfields and parks #

Any geometry where a subject cannot plausibly be is a constraint. Water is the one people notice because it is drawn in a different colour. A case in the middle of a runway or a reservoir is just as wrong and usually goes unremarked.

Edge cases or notes #

• Prepare the geometry. prep() before a loop of contains calls.
• within and contains disagree on boundaries. Points exactly on the line are a tiny but real class.
• Union the land once. Rebuilding it inside the loop dominates the runtime.
• Islands need the full multipolygon. Taking the largest part strands everyone on a small island.
• Rivers and estuaries are fiddly. A tidal river is water at some states of tide and not others.
• A constrained mask is not the mask you documented. Publish the constraint along with the radii.
• Rejection changes nothing about repeat releases. Still store the output rather than redrawing.
• Check the bounding box after masking. It grows by r_max in every direction, which can cross a coastline that was not in the input extent.

Internal links #

• How to apply donut geomasking to sensitive points in Python — the mask this fixes
• Geomasking methods explained: donut, random and adaptive — choosing a radius the geography can accommodate
• Points plot in the ocean at latitude-longitude zero — the other reason points end up in water
• How to fix a CRS mismatch in GeoPandas — the inland-water case
• How to select features by location in GeoPandas — the containment test
• Spatial k-anonymity explained — why the adaptive radius avoids the problem
• How to extract a coastline from a raster in Python — building a land mask when no vector exists
• How to run a privacy check before publishing a spatial dataset — catching this before it ships

FAQ #

Why do my geomasked points end up in the sea? #

Because the displacement is unconstrained. Any point within the outer radius of water can land in it — 23.8% of the test dataset was within 300 m of the sea, and about a tenth of those failed.

Should I snap offending points to the nearest land? #

No. That puts every rejected point on the coastline, which is visible on a map and tells an attacker the true point was inland. Redraw instead.

How many redraws does rejection sampling need? #

Few. On a coastal city dataset, 63 of 3,109 points needed a redraw, three rounds resolved all of them, and the worst single point took four attempts.

My points are in water everywhere, not just at the coast. Why? #

The mask is almost certainly running in a geographic CRS, so metres are being added to degrees. Reproject before masking.

Which coastline should I use? #

One whose accuracy is better than your masking radius. Natural Earth 10m is generalised for world maps and can be well over a hundred metres out at city scale.

Does the constraint change the privacy guarantee? #

Yes, slightly: the released point is known to be on land, which narrows the candidate set. Publish the constraint along with the radii so users and reviewers can account for it.