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 — the OpenStreetMap address points of Brighton & Hove — a 50–300 m donut mask put 83 of 4,059 points (2.04%) into the water on one draw, and a 200–1000 m mask put 607 (15.0%) 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 redrew 83 of 4,059 points and needed at most 7 attempts for a single point (Example 2). Because only the offenders are redrawn, the displacement distribution barely moves.

Count, check the CRS, measure the coastal band, redraw

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, 818 of 4,059 (20.2%) lay within 300 m of the sea, and roughly a tenth of those ended up in the water on a 50–300 m draw (Example 4):

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 — on these points, a median displacement of 40.5 m at k=10 against 213.8 m for the fixed donut (how to apply a donut geomask has the code).

A land mask and reject-and-redraw

The data used on this page

The points are every OpenStreetMap node in Brighton & Hove carrying an addr:housenumber, standing in for cases on a coastline. The sea polygon is built in Example 1 from the OpenStreetMap coastline ways in a box around the city, fetched from Overpass; the regions in Example 3 are the ONS December 2024 wards of Brighton and Hove. Overpass via the maps.mail.ru mirror.

import json
import os

import geopandas as gpd
import numpy as np
import osmnx as ox
import requests
from shapely.geometry import box

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})
addresses = features[features.index.get_level_values(0) == "node"][["geometry"]].reset_index(drop=True).to_crs(27700)
xy = np.c_[addresses.geometry.x, addresses.geometry.y]
rng = np.random.default_rng(20260915)


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


BOX_WEST, BOX_SOUTH, BOX_EAST, BOX_NORTH = gpd.GeoSeries([area], crs=4326).to_crs(27700).buffer(3000).to_crs(4326).total_bounds
box_27700 = gpd.GeoSeries([box(BOX_WEST, BOX_SOUTH, BOX_EAST, BOX_NORTH)], crs=4326).to_crs(27700).iloc[0]
if not os.path.exists("coastline.json"):
    query = (f'[out:json][timeout:120];way["natural"="coastline"]({BOX_SOUTH},{BOX_WEST},{BOX_NORTH},{BOX_EAST});'
             "out geom;")
    reply = requests.get(ox.settings.overpass_url + "/interpreter", params={"data": query}, timeout=240,
                         headers={"User-Agent": ox.settings.http_user_agent})
    reply.raise_for_status()
    open("coastline.json", "w").write(reply.text)

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.rename(columns={"WD24NM": "ward"}).to_file("wards.gpkg", driver="GPKG")
print(f"{len(addresses):,} address points; {len(json.load(open('coastline.json'))['elements'])} coastline ways; "
      f"{len(wards)} wards")

4,059 address points; 22 coastline ways; 23 wards

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_geom = sea.geometry.iloc[0]
land_geom = box_27700.difference(sea_geom)

sea polygon: 105.5 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: 83, max attempts 7

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

points = addresses
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])
print(f"{len(joined):,} points masked inside {joined.ward.nunique()} wards")

4,059 points masked inside 23 wards

Example 4 — how much of the set is at risk, by radius

from shapely import distance as sdist, points as spoints

xy = np.c_[addresses.geometry.x, addresses.geometry.y]
to_sea = sdist(spoints(xy), sea_geom)
for lo, hi in ((0, 100), (50, 300), (100, 500), (200, 1000)):
    out = xy + donut_offsets(len(xy), lo, hi, np.random.default_rng(20260915))
    wet = sdist(spoints(out), sea_geom) == 0
    band = to_sea <= hi
    print(f"{lo}-{hi} m: {wet.sum()} in water ({wet.mean():.2%} of all), "
          f"{(wet & band).sum() / band.sum():.2%} of the {band.sum():,} within {hi} m of the sea")

0-100 m: 9 in water (0.22% of all), 14.06% of the 64 within 100 m of the sea
50-300 m: 83 in water (2.04% of all), 10.15% of the 818 within 300 m of the sea
100-500 m: 263 in water (6.48% of all), 14.96% of the 1,758 within 500 m of the sea
200-1000 m: 607 in water (14.95% of all), 22.38% of the 2,712 within 1000 m of the sea

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.

Why snapping is worse than redrawing

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.

Prepared geometry, boundaries and islands

• 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, 83 of 4,059 points needed a redraw, and the worst single point took seven 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.