Spatial Epi Public Health GIS

Enforcing Spatial K-Anonymity on Case-Point Releases

Publishing a case-point map can silently violate k-anonymity in low-density areas even after coordinates are jittered, because the disc that actually contains kk people may dwarf the jitter radius. This guide, part of spatial k-anonymity for health microdata, computes the smallest population-bearing disc around each point, keeps points whose disc stays tight, and snaps the rest to a coarser geography with a suppression report you can defend.

Problem Context & Constraints

The naive point release — take each case, add a random offset of a few hundred metres, plot it — feels private because the marker no longer sits on the patient’s roof. But privacy is not measured in metres of displacement; it is measured in people. A 250 m jitter around an address in a dense urban block might enclose two thousand residents, comfortably anonymous. The identical jitter around a farmhouse might enclose four. The offset is the same; the anonymity set differs by three orders of magnitude, and only the second case is a re-identification waiting to happen.

Three constraints make point releases harder than areal counts:

  • Jitter radius and anonymity set are decoupled. A fixed displacement gives a fixed geometric uncertainty but a wildly variable population uncertainty. Uniform jitter under-protects exactly where protection matters most — the sparse periphery — and over-protects the dense core where it is least needed.
  • Density varies continuously across one map. A single global radius cannot be correct everywhere. The radius that reaches kk people is a property of the location, so it must be computed per point, not assumed.
  • Points that cannot be protected must not be dropped. A case that fails the test still happened; deleting it biases the surface toward densely populated places. The correct move is to coarsen its geography — represent it as a larger unit — not to erase it.

The fix replaces a global jitter with a per-point, population-aware rule: grow a disc around each point until it holds kk people, and if that disc is unacceptably large, stop representing the case as a point and aggregate it up to a unit that clears the threshold. This is a natural companion to donut geomasking for patient address points, which controls where a displaced point lands; here we decide whether a point may be published at all.

Population-Aware Disc: Keep or Snap Two panels. In the dense-area panel, a case point sits inside a small dashed disc that already encloses k residents, so it is released as a masked point with its precision retained. In the sparse-area panel, the disc needed to enclose k residents is far larger and exceeds the maximum allowed radius, so the case is snapped to a coarser unit such as a census tract. The lesson is that point precision is kept only where the surrounding population supports k. Reaching k people around a case point — keep or snap Dense area k people within a small radius r_i Release masked point precision retained Sparse area r_i exceeds the r_max limit Snap to coarse unit e.g. census tract Point precision is kept where population supports k; sparse points snap to a coarser geography

Prerequisites

Pin the stack so the smallest-disc computation and the resulting release are reproducible:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x, pyproj 3.6.x)
  • numpy 1.26.4, scipy 1.13.1 — cKDTree for radius queries

Input state assumed below:

  • Case points as a point GeoDataFrame with stable, non-patient identifiers and a defined CRS.
  • Population basis as blocks (or a raster converted to weighted centroids) carrying the at-risk count matching the case definition — the same denominator discipline the parent section insists on.
  • CRS: one projected, area-appropriate CRS for all layers so a radius is measured in metres. Reproject before any distance query; getting WGS84 and a metric grid aligned first is covered in coordinate reference systems for public health.
  • A declared kmin⁡k_{\min} and a maximum acceptable radius rmax⁡r_{\max}, both from written policy. rmax⁡r_{\max} is the line past which a “point” is too vague to be worth publishing as a point.

Step-by-Step Solution

For each case point, sort nearby population cells by distance, accumulate their at-risk population outward, and record rir_i as the radius at which the running total first reaches kmin⁡k_{\min}. If ri≤rmax⁡r_i \le r_{\max} the point is releasable (optionally masked within rir_i); if ri>rmax⁡r_i > r_{\max} the case is snapped to the coarse unit that contains it. Every decision is logged, and a suppression report tallies what was coarsened.

# Per-point smallest-disc-to-k, with snap-to-coarse fallback and a suppression report.
# python 3.11
# geopandas==0.14.4  shapely==2.0.4  pyproj==3.6.1
# numpy==1.26.4      scipy==1.13.1
import logging
import numpy as np
import geopandas as gpd
from scipy.spatial import cKDTree

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("k_point")

EQUAL_AREA = "EPSG:5070"  # metric, area-faithful — radii are in metres

def smallest_disc_radius(case_xy, pop_tree, pop_xy, pop_at_risk, k_min, r_hard=50_000):
    """
    Radius at which cumulative at-risk population around a case first reaches k_min.
    Queries neighbors out to a hard cap; returns np.inf if k_min is unreachable.
    """
    # candidate cells within the hard cap, nearest first
    idx = pop_tree.query_ball_point(case_xy, r=r_hard)
    if not idx:
        return np.inf
    d = np.linalg.norm(pop_xy[idx] - case_xy, axis=1)
    order = np.argsort(d, kind="stable")           # deterministic ordering
    cum = np.cumsum(pop_at_risk[idx][order])
    reach = np.searchsorted(cum, k_min)             # first index where cum >= k_min
    if reach >= len(cum):
        return np.inf                               # not enough people even at the cap
    return float(d[order][reach])

def enforce_point_k(cases, pop, coarse, k_min, r_max, geoid="geoid",
                    pop_field="pop_at_risk", coarse_id="tract_id"):
    """
    cases : point GeoDataFrame (release candidates)
    pop   : population basis (points or block centroids) with pop_field
    coarse: coarse polygons (e.g. tracts) with coarse_id — the snap target
    Returns (released_points, aggregated_report).
    """
    for name, gdf in (("cases", cases), ("pop", pop), ("coarse", coarse)):
        if gdf.crs is None:
            raise ValueError(f"{name} has no CRS; refuse to measure radii on undefined units.")
    cases = cases.to_crs(EQUAL_AREA).sort_values(geoid).reset_index(drop=True)
    pop = pop.to_crs(EQUAL_AREA)
    coarse = coarse.to_crs(EQUAL_AREA)

    pop_xy = np.column_stack([pop.geometry.x, pop.geometry.y])
    pop_pts = pop[pop_field].to_numpy(float)
    tree = cKDTree(pop_xy)

    radii, keep_mask = [], []
    for _, row in cases.iterrows():
        cxy = np.array([row.geometry.x, row.geometry.y])
        r_i = smallest_disc_radius(cxy, tree, pop_xy, pop_pts, k_min)
        radii.append(r_i)
        keep_mask.append(r_i <= r_max)              # releasable as a point?
    cases["disc_radius_m"] = radii
    cases["releasable"] = keep_mask

    released = cases[cases["releasable"]].copy()
    to_snap = cases[~cases["releasable"]].copy()

    # snap unreachable cases to the coarse unit that contains them, then count per unit
    if len(to_snap):
        snapped = gpd.sjoin(to_snap, coarse[[coarse_id, "geometry"]],
                            predicate="within", how="left")
        report = (snapped.groupby(coarse_id)
                         .size().rename("aggregated_cases").reset_index())
    else:
        report = coarse[[coarse_id]].iloc[0:0].assign(aggregated_cases=0)

    log.info("released_points=%d snapped_cases=%d coarse_units=%d k_min=%d r_max=%.0fm",
             len(released), len(to_snap), report.shape[0], k_min, r_max)
    return released, report

The contract mirrors the areal case: nothing is deleted and nothing is imputed. A point that clears rmax⁡r_{\max} is published (optionally after masking anywhere inside its own rir_i, which is guaranteed anonymous by construction); a point that fails is represented honestly as a coarse-unit count. The disc_radius_m column is retained because it is the achieved-privacy measurement — a reviewer can confirm the guarantee from the output alone.

Validation & Edge Cases

A released layer needs three assertions before it leaves the pipeline.

1. Minimum-k assertion on released points. By construction every released point had ri≤rmax⁡r_i \le r_{\max}, but a bug in the population join can quietly zero the denominator. Re-derive the enclosed population for a sample and assert it clears the threshold:

def assert_min_k(released, pop, k_min, pop_field="pop_at_risk"):
    tree = cKDTree(np.column_stack([pop.geometry.x, pop.geometry.y]))
    pop_pts = pop[pop_field].to_numpy(float)
    worst = np.inf
    for _, row in released.iterrows():
        idx = tree.query_ball_point([row.geometry.x, row.geometry.y], r=row["disc_radius_m"] + 1)
        worst = min(worst, float(pop_pts[idx].sum()))
    assert worst >= k_min, f"released point below k_min: enclosed pop {worst:.0f} < {k_min}"
    return worst

2. Count of aggregated points reconciles. Released points plus snapped cases must equal the input count exactly — a mismatch means a case fell through a null join and vanished:

INFO enforce_point_k: released_points=4127 snapped_cases=619 coarse_units=88 k_min=5 r_max=2000m
assert len(released) + int(report["aggregated_cases"].sum()) == len(cases), \
    "case accounting mismatch: a point was dropped, not coarsened"

3. Snap target actually contains the point. A left join to coarse polygons produces null tract_id for any case that falls in a gap between polygons (a coastline sliver, a boundary crack). Those nulls must be caught and routed to review, never released as an unlabeled aggregate. A wall of nulls almost always signals a CRS mismatch or an incomplete coarse layer rather than genuinely off-map cases.

Two further failure modes recur in practice. First, an unreachable kmin⁡k_{\min} returns ri=∞r_i = \infty — the case sits in a region so sparse that even the hard query cap holds fewer than kk people; these must snap to coarse geography and, if the coarse unit itself is below threshold, escalate to the areal aggregation in the parent guide. Second, clustered cases inflate false confidence: two cases at the same address each individually reach kk, but publishing both points at once narrows the joint anonymity set, so co-located cases should share one masked location or be aggregated together.

Compliance Notes

Three artifacts make a point release defensible:

  • The per-point radius, retained in the output. disc_radius_m is the auditable proof that each published point met the population threshold; it lets a reviewer reconstruct the guarantee without re-running the whole pipeline, and because it varies per point it also documents where the map is precise and where it is deliberately vague.
  • The suppression report. The count of cases snapped to each coarse unit is the record of what was not published as a point. That tally is the difference between an honest release and one that quietly biases toward dense areas by deleting sparse cases.
  • Pinned parameters and vintage. Store kmin⁡k_{\min}, rmax⁡r_{\max}, the population layer identity and vintage, and library versions alongside the release. Two runs under different rmax⁡r_{\max} values are different disclosure decisions and must not be compared as if equivalent.

Once a masked point layer exists, close the loop by re-measuring it: the achieved anonymity of the final map should be confirmed with the same population-based test used across this section, as detailed in validating geomasks with spatial k-anonymity. A mask you have not validated against a denominator is an assumption, not a guarantee.