How to add differentially private noise to spatial counts

Problem statement

Suppression answers "is this table safe?" once. If you publish the same geography every month for five years, you need an answer that composes — one where running the query again has a cost you can see. That is what a differential privacy budget gives you, and the Laplace mechanism is the three-line implementation.

The code is easy and the accounting is not. Three things have to be decided before any noise is added: what one unit of privacy is, how many cells one unit can affect, and how much total budget the product gets over its lifetime. This guide implements the mechanism on 51 cells of 500 m London crime counts and shows what each choice does to the map.

Quick answer

import numpy as np

def private_counts(counts, epsilon, sensitivity=1, rng=None):
    """Laplace mechanism: noise scale = sensitivity / epsilon."""
    rng = rng or np.random.default_rng()
    noisy = counts + rng.laplace(0.0, sensitivity / epsilon, size=len(counts))
    return np.maximum(np.round(noisy), 0).astype(int)      # post-processing is free

released = private_counts(cell_counts.values, epsilon=1.0, rng=np.random.default_rng(7))
print(f"mean absolute error {np.abs(released - cell_counts.values).mean():.2f} against a median count of "
      f"{cell_counts.median():.0f}")

mean absolute error 1.10 against a median count of 211

Measured on the real counts, ε=1.0 gave a mean absolute error of 1.10 crimes per cell against a median cell count of 211 — invisible in any use of the map. At ε=0.1 the mean error in Example 1 was 10.81 and the 95th percentile 41.83; no cell went negative only because the emptiest cell here held 33 crimes.

Step-by-step solution

1. State the unit of privacy

Write it in the metadata in words: "adding or removing all records of one person changes the output distribution by at most $e^\varepsilon$". If you cannot say that sentence truthfully — because one person contributes many records and you have not bounded it — the ε you publish is not the ε you have.

2. Bound the sensitivity

Sensitivity is the maximum change one unit can cause across all published cells. For counts where one person appears once, it is 1. For repeat visits or trajectories, clamp the contribution:

def clamp(df, subject, k):
    return df.groupby(subject, group_keys=False).head(k)

clamped = clamp(visits, "patient_id", k=3)     # sensitivity is now 3

The clamp itself is a modelling choice with a bias — it truncates the heaviest users — and belongs in the documentation.

3. Choose ε and write down the budget

ε is a policy parameter. Decide the total for the product's lifetime, then divide it among the queries you intend to publish. A ledger of two columns — query, ε spent — is enough, and it must include the queries you ran to design the product.

4. Add the noise to every cell, including the empty ones

A cell that is omitted because its true count is zero leaks that the count is zero. Build the full cell list from the geography, reindex the counts onto it with zeros, and noise all of them.

counts = observed.reindex(all_cells, fill_value=0)
noisy = private_counts(counts.values, epsilon=1.0)

5. Post-process for usability

Clamping negatives, rounding to integers, smoothing and reconciling a hierarchy so children sum to a parent are all allowed, because none of them looks at the raw data again.

6. Report the accuracy you delivered

Users need to know the noise scale to interpret the map. Publish ε, the sensitivity, the noise distribution and the resulting error percentiles at the cell sizes you offer.

7. Check whether the map survived

If the interesting features are small counts, the mechanism will destroy them and no ε you can defend will save them. At that point, aggregate to bigger units instead — which is a different, non-private-but-adequate answer, and an honest one.

Code examples

The data used on this page

Street-level crime from data.police.uk for January to March 2025, for a box over Westminster, Holborn and Bloomsbury (0.16°W to 0.10°W, 51.505°N to 51.530°N), fetched in tiles because the API returns at most 10,000 crimes a request. The cells are 500 m squares aligned to the British National Grid, and the ward counts are the same crimes joined to the ONS December 2024 wards of the City of Westminster. police.uk moves every crime to one of a set of anonymised map points, so the counts are counts at those points.

import os

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

WEST, SOUTH, EAST, NORTH = -0.16, 51.505, -0.10, 51.530      # Westminster, Holborn, Bloomsbury
if not os.path.exists("london_crime_2025q1.csv"):
    rows = []
    for month in ("2025-01", "2025-02", "2025-03"):
        for x0, x1 in zip(np.linspace(WEST, EAST, 4)[:-1], np.linspace(WEST, EAST, 4)[1:]):
            for y0, y1 in zip(np.linspace(SOUTH, NORTH, 3)[:-1], np.linspace(SOUTH, NORTH, 3)[1:]):
                reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                      data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": month})
                reply.raise_for_status()                    # a 503 means over 10,000 crimes: split further
                rows += [{"id": c["id"], "month": c["month"], "category": c["category"],
                          "lon": float(c["location"]["longitude"]), "lat": float(c["location"]["latitude"]),
                          "street_id": c["location"]["street"]["id"]} for c in reply.json()]
    pd.DataFrame(rows).drop_duplicates("id").to_csv("london_crime_2025q1.csv", index=False)
crimes = pd.read_csv("london_crime_2025q1.csv")
crimes = crimes[crimes.lon.between(WEST, EAST) & crimes.lat.between(SOUTH, NORTH)]
crimes = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
print(f"{len(crimes):,} crimes, {crimes.category.nunique()} categories, {crimes.street_id.nunique():,} snap points")

x, y = crimes.geometry.x, crimes.geometry.y
cell_counts = pd.Series(list(zip((x // 500).astype(int), (y // 500).astype(int)))).value_counts().sort_index()

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query")
wards_gdf = gpd.read_file(requests.get(ONS, params={"where": "LAD24NM='Westminster'", "outFields": "WD24CD,WD24NM",
                                                    "outSR": 27700, "f": "geojson"}, timeout=120).text)
in_ward = gpd.sjoin(crimes, wards_gdf, predicate="within")
ward_counts = in_ward.groupby("WD24NM").size().reindex(wards_gdf.WD24NM, fill_value=0).to_numpy()
print(f"{len(cell_counts)} cells of 500 m: {cell_counts.sum():,} crimes, median {cell_counts.median():.0f}, "
      f"min {cell_counts.min()}, max {cell_counts.max():,}; {len(wards_gdf)} Westminster wards, {ward_counts.sum():,} crimes")

22,293 crimes, 14 categories, 1,249 snap points
51 cells of 500 m: 22,293 crimes, median 211, min 33, max 3,533; 18 Westminster wards, 15,686 crimes

Example 1 — the accuracy/ε curve on real counts

import numpy as np

rng = np.random.default_rng(7)
counts = cell_counts.values          # 500 m cells, central London, 2025 Q1

print(f"{'eps':>5} {'mean |err|':>10} {'p95 |err|':>9} {'negatives':>10} {'true>0 -> 0':>12} {'drift':>8}")
for eps in (0.1, 0.5, 1.0, 2.0):
    noisy = counts + rng.laplace(0, 1 / eps, len(counts))
    err = np.abs(noisy - counts)
    clipped = np.maximum(np.round(noisy), 0)
    print(f"{eps:5.1f} {err.mean():10.2f} {np.percentile(err, 95):9.2f} "
          f"{int((noisy < 0).sum()):10,} {int(((counts > 0) & (clipped == 0)).sum()):12,} "
          f"{clipped.sum() - counts.sum():+8.0f}")

  eps mean |err| p95 |err|  negatives  true>0 -> 0    drift
  0.1      10.81     41.83          0            0      -53
  0.5       2.00      5.55          0            0       -4
  1.0       0.94      2.79          0            0       +9
  2.0       0.53      1.42          0            0       -5

The smallest true count in these cells was 33 and the largest 3,533, so ε=1 leaves the map visually identical and ε=0.1 does not — and a cell of 2 or 3, common at finer resolutions, would be noise at ε=0.1.

Example 2 — a budget ledger you can actually keep

import json, pathlib, datetime

class Budget:
    def __init__(self, total, path="epsilon_ledger.json"):
        self.total, self.path = total, pathlib.Path(path)
        self.spent = json.loads(self.path.read_text()) if self.path.exists() else []

    def remaining(self):
        return self.total - sum(e["epsilon"] for e in self.spent)

    def spend(self, label, epsilon):
        if epsilon > self.remaining():
            raise ValueError(f"{label}: {epsilon} requested, {self.remaining():.3f} left")
        self.spent.append({"label": label, "epsilon": epsilon,
                           "when": datetime.datetime.now().isoformat(timespec="seconds")})
        self.path.write_text(json.dumps(self.spent, indent=2))
        return epsilon

pathlib.Path("epsilon_ledger.json").unlink(missing_ok=True)   # start the example with an empty ledger
budget = Budget(total=4.0)
eps = budget.spend("2025-Q1 500 m crime map", 1.0)
released = private_counts(counts, eps)
print(f"remaining budget: {budget.remaining():.2f}")

remaining budget: 3.00

The ledger is the artefact that makes the guarantee real. Without it, "differentially private" describes one run, not a product.

Example 3 — reconciling a hierarchy after noising

import numpy as np

def reconcile(child_counts, parent_total):
    """Scale noisy children so they sum to a noisy parent, then round to integers."""
    child = np.maximum(child_counts, 0).astype(float)
    if child.sum() == 0:
        return np.zeros_like(child, dtype=int)
    scaled = child * (parent_total / child.sum())
    out = np.floor(scaled).astype(int)
    short = int(parent_total - out.sum())
    if short > 0:                                   # hand the remainder to the largest fractions
        order = np.argsort(-(scaled - out))[:short]
        out[order] += 1
    return out

wards = private_counts(ward_counts, epsilon=0.5)
borough = private_counts(np.array([ward_counts.sum()]), epsilon=0.5)[0]
print(reconcile(wards, borough).sum(), borough)

15684 15684

Both queries spend budget; the reconciliation itself does not, because it only touches the two noisy outputs.

Explanation

Why the noise does not depend on the count

The Laplace scale is sensitivity/ε, and sensitivity for a count is how much one person can change it — one. That is independent of whether the cell holds 3 or 3,533. It makes the mechanism nearly free on large counts and fatal on small ones, which is the single most important practical fact about it.

Why every cell must be published

If cells appear in the output only when their true count is non-zero, the presence of a cell is itself a noiseless statistic. Build the cell list from the geography, not from the data.

Why post-processing is free

The guarantee is a property of the mapping from the raw data to the released values. Anything computed from the released values alone cannot increase the information about the raw data, so clamping, rounding and smoothing are safe — and they should be used, because raw Laplace output is unusable.

Why this is often the wrong tool for a map

Maps are interesting where counts are small: the one cluster, the rare category, the quiet ward. Those are exactly the cells noise destroys. If your map's purpose survives ε=1 on 500 m cells, use it; if the purpose is to find sparse events, aggregation to units built for the threshold is the better answer. How to aggregate points to units that meet a minimum count covers that route.

Edge cases or notes

• Do not resample negative values. Clamp them; resampling is data-dependent and breaks the guarantee.
• Do not suppress noisy cells. Suppression after noising reintroduces a raw-data decision.
• The geometry is public. Only the counts are protected.
• One person, many cells. Trajectories need a clamp on cells touched, not just on records.
• Seeds are for reproducibility, not secrecy. Publishing the seed lets a user reproduce the file; it does not weaken ε.
• Composition is additive for pure ε. Advanced composition gives a better bound but needs (ε, δ).
• Budget the design queries. Choosing a cell size by looking at counts costs ε.
• Small ε makes rates nonsense. A noisy numerator over a small denominator can exceed 100%.

Internal links

• Differential privacy for spatial counts explained — what ε means and where the mechanism comes from
• How to suppress small counts in a spatial table in Python — the alternative, and why the two do not mix
• How to aggregate points to units that meet a minimum count — when the map needs the sparse cells
• Aggregation and suppression rules explained — thresholds and margins
• Counts, rates and densities explained — what noise does to a denominator
• How to smooth rates with empirical Bayes in Python — post-processing noisy small cells
• How to run a privacy check before publishing a spatial dataset — where the ledger check belongs
• Geoprivacy explained: why coordinates are personal data — choosing between the defences

FAQ

How much noise should I add?

Scale sensitivity/ε. For a count where one person appears once and ε=1, that is Laplace with scale 1 — about one count of error per cell.

What ε should I use?

There is no correct value; it is a policy choice. Published products range from about 0.1 to 10 in total. Decide the lifetime budget first, then divide it.

Do I add noise to zero-count cells?

Yes. Omitting them publishes the fact that they are zero without any noise at all.

Can I clamp negative counts to zero?

Yes. Post-processing a released value never weakens the guarantee, and raw Laplace output is unusable without it.

What if one person contributes many records?

Clamp to at most k records per person and set sensitivity to k. Without the clamp, one heavy user sets the noise scale for everyone.

Can I combine differential privacy with suppression?

No. Suppressing cells based on their values is a raw-data decision that is not covered by the budget.