How to measure re-identification risk in a point dataset

Problem statement

Before you can decide how to protect a point dataset you have to know how exposed it is, and "it looks fine on the map" is not a measurement. The audit has to produce numbers somebody can put in a data-sharing agreement: how many records sit alone on a location, how many candidate subjects each released point leaves, what the smallest publishable cell contains, and which handful of records carry most of the risk.

The hard part is not the code. It is remembering that the dataset you are auditing is not the file — it is the file plus every filter the release exposes, plus whatever an attacker can join to it. This guide runs that audit end to end on real data: 4,064 OpenStreetMap address points in Brighton & Hove, and 4,013 UK police crime records for central Leicester that have already been anonymised by the publisher, audited against Ordnance Survey's complete list of 39,262 addressable locations in the same box.

Quick answer

Three numbers, in this order:

import numpy as np, pandas as pd, geopandas as gpd
from scipy.spatial import cKDTree

pts = gpd.read_file("release.gpkg").to_crs(27700)          # a metric CRS
xy = np.c_[pts.geometry.x, pts.geometry.y]

# 1. how unique is each location?
key = pd.Series(list(map(tuple, xy.round(1))))
print("rows alone on their coordinate:", int((key.value_counts() == 1).sum()))

# 2. how many candidate subjects does each point leave?
tree = cKDTree(np.c_[addresses.geometry.x, addresses.geometry.y])
k = np.array([len(tree.query_ball_point(p, 50)) for p in xy])
print(f"k at 50 m: min {k.min()}, median {np.median(k):.0f}, below 5: {(k < 5).sum()}")

# 3. what is the smallest cell the release exposes?
cells = pts.groupby([key, pts["month"], pts["category"]]).size()
print("smallest published cell:", cells.min())

rows alone on their coordinate: 51
k at 50 m: min 1, median 68, below 5: 222
smallest published cell: 1

On the published Leicester crime extract, already snapped to 368 shared points by the police, only 51 records are alone on their coordinate and the median point has 68 addresses within 50 m — but 222 points have fewer than five, and the smallest published cell is 1, in 1,501 of 2,211 cells. On the unmasked Brighton addresses the first number is 4,064 of 4,064: every record alone.

Metres first, then uniqueness, candidates and filters

1. Reproject to metres before anything else

Every distance in the audit is in metres. Reproject once at the top, and record the CRS in the report — a risk number without a CRS is not reproducible.

2. Measure coordinate uniqueness

Count how many records share each rounded coordinate. Round at the precision you intend to publish, not at full precision, or you measure the source rather than the release.

for dp in (6, 5, 4, 3, 2):
    key = pd.Series(list(zip(lat.round(dp), lon.round(dp))))
    vc = key.value_counts()
    print(f"{dp} dp: {len(vc):5,} distinct, {int((vc == 1).sum()):5,} occupied by one record")

6 dp: 4,064 distinct, 4,064 occupied by one record
5 dp: 4,063 distinct, 4,062 occupied by one record
4 dp: 3,285 distinct, 2,603 occupied by one record
3 dp:   754 distinct,   246 occupied by one record
2 dp:    72 distinct,     4 occupied by one record

Three decimal places — a cell of about 111 m by 70 m at this latitude — still leaves 246 of the 4,064 points alone in their cell.

3. Count candidate subjects against a real auxiliary file

This is the number that means something. Use the most complete address or premises file you can obtain, and record which one. Spatial k-anonymity explained covers what the count means.

4. Enumerate the smallest cell the release exposes

Cross every dimension a user can filter by. The Leicester extract is a good example of the gap between the two:

5. Rank the records by risk and look at the top

Risk is concentrated. The isolated points, the rare categories and the records whose attributes are unique in their neighbourhood are a small, enumerable set, and they are where the release decision actually lives.

6. Test the attributes as well as the geometry

For each record, count how many other records share both its neighbourhood and its attribute combination. A record that is unique on both is a re-identification waiting for an auxiliary file.

7. Write the report as a file, not a message

The audit belongs next to the data, with the date, the auxiliary source, the CRS, the parameters and the numbers. How to run a privacy check before publishing a spatial dataset turns it into a pipeline gate.

A uniqueness and candidate-count audit

The data used on this page

Three real extracts. Every OpenStreetMap node in Brighton & Hove carrying an addr:housenumber (the uniqueness profile in step 2, as lat/lon); street-level crime from data.police.uk for January to March 2025 in a box over central Leicester (1.150°W to 1.115°W, 52.625°N to 52.645°N), saved as leicester_2025Q1.gpkg and as the release.gpkg of the quick answer; and, as the attacker's auxiliary file (addresses), every UPRN in the same Leicester box from Ordnance Survey's OS Open UPRN (the September 2026 release, 619 MB zipped, Open Government Licence) — the complete, free list of addressable locations in Great Britain. OpenStreetMap holds only 431 address nodes in that box, which would make almost every point look like it had no candidates at all. Overpass via the maps.mail.ru mirror, read from its answer of 9 October 2026, which the geoprivacy guides share. The crimes are a copy of the API's answer saved on 11 October 2026.

import io
import os
import time
import zipfile

import geopandas as gpd
import numpy as np
import osmnx as ox
import pandas as pd
import requests

# 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-measure-reidentification-risk-python-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})
brighton = features[features.index.get_level_values(0) == "node"].reset_index(drop=True)
lat, lon = brighton.geometry.y, brighton.geometry.x

WEST, SOUTH, EAST, NORTH = -1.150, 52.625, -1.115, 52.645        # central Leicester
if not os.path.exists("leicester_2025Q1.gpkg"):
    rows = []
    # live: requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120, data={"date": month,
    #     "poly": f"{SOUTH},{WEST}:{SOUTH},{EAST}:{NORTH},{EAST}:{NORTH},{WEST}"}) for each month, retrying on 429;
    # read instead as the API answered on 11 October 2026, because police.uk revises past months
    saved = requests.get("https://data.spatialworkflow.io/police-leicester-2025q1-2026-10-11.json", timeout=120).json()
    for month in ("2025-01", "2025-02", "2025-03"):
        rows += [{"id": c["id"], "month": c["month"], "category": c["category"],
                  "lon": float(c["location"]["longitude"]), "lat": float(c["location"]["latitude"])}
                 for c in saved[month]]
    # sorted by id: the API returns the same crimes in a different order on each call
    table = pd.DataFrame(rows).drop_duplicates("id").sort_values("id")
    gpd.GeoDataFrame(table, geometry=gpd.points_from_xy(table.lon, table.lat), crs=4326).to_file("leicester_2025Q1.gpkg")
gpd.read_file("leicester_2025Q1.gpkg").to_file("release.gpkg")

if not os.path.exists("leicester_uprn.gpkg"):
    if not os.path.exists("osopenuprn_csv.zip"):
        with requests.get("https://api.os.uk/downloads/v1/products/OpenUPRN/downloads"
                          "?area=GB&format=CSV&redirect", stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open("osopenuprn_csv.zip", "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)
    import zipfile
    with zipfile.ZipFile("osopenuprn_csv.zip") as z:
        name = next(n for n in z.namelist() if n.endswith(".csv"))
        with z.open(name) as f:
            keep = [c[c.LONGITUDE.between(WEST, EAST) & c.LATITUDE.between(SOUTH, NORTH)]
                    for c in pd.read_csv(f, chunksize=2_000_000)]
    uprn = pd.concat(keep)
    gpd.GeoDataFrame(uprn[["UPRN"]], crs=27700,
                     geometry=gpd.points_from_xy(uprn.X_COORDINATE, uprn.Y_COORDINATE)).to_file("leicester_uprn.gpkg")
addresses = gpd.read_file("leicester_uprn.gpkg")
print(f"Brighton & Hove: {len(brighton):,} address nodes; Leicester box: {len(addresses):,} UPRNs")

Brighton & Hove: 4,064 address nodes; Leicester box: 39,262 UPRNs

Example 1 — the full uniqueness profile

import numpy as np, pandas as pd, geopandas as gpd

def uniqueness_profile(gdf, quasi_identifiers=(), round_m=1):
    xy = pd.Series(list(zip((gdf.geometry.x / round_m).round().astype(int),
                            (gdf.geometry.y / round_m).round().astype(int))))
    out = {"rows": len(gdf), "distinct_locations": xy.nunique()}
    vc = xy.value_counts()
    out["locations_with_one_record"] = int((vc == 1).sum())
    if quasi_identifiers:
        key = pd.Series(list(zip(xy, *[gdf[c] for c in quasi_identifiers])))
        kvc = key.value_counts()
        out["quasi_id_cells"] = len(kvc)
        out["cells_of_one"] = int((kvc == 1).sum())
        out["records_in_cells_of_one"] = int(kvc[kvc == 1].sum())
    return out

crimes = gpd.read_file("leicester_2025Q1.gpkg").to_crs(27700)
print(uniqueness_profile(crimes, ["month", "category"]))

{'rows': 4013, 'distinct_locations': 368, 'locations_with_one_record': 51, 'quasi_id_cells': 2211, 'cells_of_one': 1501, 'records_in_cells_of_one': 1501}

Example 2 — candidate counts at several radii

from scipy.spatial import cKDTree

tree = cKDTree(np.c_[addresses.geometry.x, addresses.geometry.y])
released = np.c_[pts.geometry.x, pts.geometry.y]

rows = []
for radius in (25, 50, 100, 250, 500):
    k = np.array([len(tree.query_ball_point(p, radius)) for p in released])
    rows.append({"radius_m": radius, "k_min": int(k.min()),
                 "k_median": float(np.median(k)),
                 "below_5": int((k < 5).sum()), "below_10": int((k < 10).sum())})
print(pd.DataFrame(rows).to_string(index=False))

 radius_m  k_min  k_median  below_5  below_10
       25      0      14.0      888      1396
       50      1      68.0      222       355
      100      4     310.0        3         3
      250    113    1889.0        0         0
      500   1265    7298.0        0         0

Report the whole table. A single radius hides the fact that the failing points are the same ones at every radius.

Example 3 — the records that carry the risk

d, _ = cKDTree(released).query(released, k=6)
isolation = d[:, 5]                                  # distance to the 5th nearest released point

worst = pts.assign(isolation_m=isolation).nlargest(20, "isolation_m")
print(worst[["category", "isolation_m"]].round(0).to_string(index=False))
print(f"\n{(isolation > 150).sum()} points are more than 150 m from their 5th nearest neighbour")

             category  isolation_m
        violent-crime        228.0
         public-order        218.0
        violent-crime        186.0
        violent-crime        179.0
        violent-crime        179.0
          other-theft        179.0
         public-order        170.0
        violent-crime        152.0
        vehicle-crime        152.0
                drugs        152.0
         public-order        141.0
        violent-crime        136.0
anti-social-behaviour        136.0
          other-crime        136.0
        violent-crime        136.0
anti-social-behaviour        136.0
          other-crime        136.0
        violent-crime        135.0
         public-order        135.0
        violent-crime        133.0

10 points are more than 150 m from their 5th nearest neighbour

Example 5 — the same audit against a different auxiliary file

features = ox.features_from_bbox((WEST, SOUTH, EAST, NORTH), tags={"addr:housenumber": True})
osm = features[features.index.get_level_values(0) == "node"].to_crs(27700)
for name, aux in [("OS Open UPRN", addresses), ("OpenStreetMap", osm)]:
    k = np.array([len(c) for c in cKDTree(np.c_[aux.geometry.x, aux.geometry.y]).query_ball_point(released, 50)])
    print(f"{name:13s} {len(aux):6,} addresses; crimes with fewer than 5 within 50 m: {(k < 5).sum():,}")

OS Open UPRN  39,262 addresses; crimes with fewer than 5 within 50 m: 222
OpenStreetMap    431 addresses; crimes with fewer than 5 within 50 m: 3,228

Example 4 — the unit table in step 4

snap = pd.Series(list(zip(crimes.geometry.x.round(1), crimes.geometry.y.round(1))), index=crimes.index)
for name, cols in [("snap point", []), ("point × month", ["month"]),
                   ("point × category", ["category"]), ("point × month × category", ["month", "category"])]:
    cells = crimes.groupby([snap, *[crimes[c] for c in cols]]).size()
    ones = int((cells == 1).sum())
    print(f"{name:26s} {len(cells):6,} {ones:6,} ({ones / len(cells):.1%})")

snap point                    368     51 (13.9%)
point × month                 866    259 (29.9%)
point × category            1,517    840 (55.4%)
point × month × category    2,211  1,501 (67.9%)

Twenty rows is a list a person can read and make a decision about — drop them, generalise them, or accept them. That is what the audit is for.

Measure at the published precision, against a named file

Why uniqueness is measured at the published precision

The source file is always unique; that is not news. What matters is whether the release is unique, so every uniqueness count must be computed after rounding, masking or aggregation, using exactly the values you intend to write.

Why the auxiliary file must be named

The candidate count is a ratio between your points and somebody else's records. In the central Leicester box OpenStreetMap holds 431 address nodes and OS Open UPRN 39,262: against OpenStreetMap, 3,228 of the 4,013 crimes had fewer than five candidates within 50 m; against UPRN, 222 did. The choice of file moved the answer by an order of magnitude, so it has to be stated.

Why the filter cross-product is the real unit

Users do not download the file; they filter the map. The Leicester extract is published as an already-anonymised product and is safe at the level of a bare point. At the level its own map exposes, two-thirds of cells describe one incident, and a person who knows roughly when and what will recover it.

Why risk is reported as a distribution

A mean candidate count of 400 tells you nothing about the point with two candidates, and the point with two candidates is the release. Always report the minimum, the count below the threshold, and the list of offenders.

Audit outputs, empty attributes and duplicates

• Audit the output of your pipeline, not the input. Masking code with a bug produces an unmasked release that the input-side audit passed.
• Empty attributes are attributes. A null category can be as distinguishing as a rare one.
• Duplicated subjects inflate k. Two records for the same person are one candidate.
• Check the CRS on both layers. A silent CRS mismatch makes every distance meaningless and every k enormous.
• Sample only for speed, never for the report. The rare points are the ones a sample misses.
• Re-run after every schema change. Adding a column changes the quasi-identifier set.
• Keep the audit script in the repository. It is the evidence, and it has to be re-runnable.
• Include the derived products. A KDE raster and an SVG map are part of the release.

Internal links

• Re-identification risk in spatial data explained — what the numbers mean
• Spatial k-anonymity explained — the candidate-count threshold
• How to run a privacy check before publishing a spatial dataset — turning the audit into a gate
• Geoprivacy explained: why coordinates are personal data — why this audit exists
• Anonymised points still identify individuals — what a passing audit still misses
• How to find the nearest point in GeoPandas — the neighbour search underneath
• How to validate a GeoDataFrame against a schema before analysis — the same idea for data quality
• How to reduce coordinate precision safely in Python — acting on the uniqueness profile

FAQ

What should I measure first?

How many records are alone on their coordinate at the precision you intend to publish. If that number is large, nothing else in the audit matters yet.

Which address file should I use as the attacker's data?

The most complete one you can legally obtain, because that is what an attacker will use. In Great Britain that is OS Open UPRN, which is free. OpenStreetMap is patchy: in central Leicester it held about one address node for every 90 UPRNs, and it made the release look far riskier than it is. That is conservative, but it is not a measurement.

What radius should I count candidates in?

The uncertainty your release admits to: the rounding cell, the masking radius, or the aggregation polygon. Report several radii so the shape of the risk is visible.

Why does my already-anonymised dataset still fail?

Usually because the threshold was applied to the map and not to the filters. Crossing a safe-looking snap point with month and category made 67.9% of cells singletons in the Leicester extract.

Do attributes need auditing too?

Yes. Count how many records share both a neighbourhood and an attribute combination; uniqueness on both is the condition an attacker needs.

How often should the audit run?

On every release, automatically, and again whenever a column is added or a filter is exposed in the UI.