Differentially Private Spatial Aggregation

A city transportation agency wants to publish a live density heatmap of shared-scooter pickups so that an analytics dashboard can show planners where demand concentrates by hour — but the raw pickup log is a book of home addresses, first-trip origins, and 6 a.m. departure points that resolve to single dwellings. Differentially private spatial aggregation is the release discipline that makes that publication safe: bin the points into a fixed grid of cells, add calibrated noise to every cell’s count, and repair the two properties the noise breaks — negative counts and a total that no longer matches. Done correctly, the dashboard sees a density surface whose per-cell error is bounded and whose privacy guarantee holds under a declared $\varepsilon$; done carelessly, the same pipeline manufactures phantom hotspots in empty parkland and leaks the exact tile where a rider lives.

This guide sits inside the Differential Privacy for Geospatial Data section and focuses on the aggregate release path — counts and densities over a partitioned plane — rather than perturbing individual coordinates. The noise itself is drawn from the mechanisms catalogued in spatial noise mechanisms; here the concern is everything around the noise: how the choice of grid resolution trades off against $\varepsilon$, why sparse regions are where the release fails first, and how uniform binning differs from the data-adaptive decompositions covered in private quadtree decomposition for spatial histograms. When the end product is a rendered surface rather than a table of cells, continue to building differentially private heatmaps, which layers interpolation and colour-scaling on top of the grid produced here.

Key concept. A uniform grid is a data-independent partition, so choosing it costs nothing from the privacy budget — every $\varepsilon$ you spend goes into the counts. That is its great advantage over adaptive decompositions, and its great weakness: the same cell size that keeps dense downtown bins accurate drowns sparse suburban bins in noise. Grid resolution is therefore the first and most consequential privacy parameter you set.

Prerequisites & Accounting Environment

The reference implementation depends only on numpy for the noise and histogram arithmetic and geopandas/shapely for reading points and emitting cell polygons. No cryptographic stack is required — unlike secret sharing or homomorphic encryption, aggregate DP runs entirely inside a single trusted curator that already holds the raw points and is trusted to add the noise honestly. If your deployment cannot make that trust assumption, you are in the local-DP regime instead, and the per-user perturbation there is a different pipeline entirely.

Three invariants must hold before the first release:

• Bounded per-user contribution. The $L_1$ sensitivity of a histogram is the maximum number of cells one data subject can change, times the magnitude of that change. If every rider contributes exactly one pickup to exactly one cell, sensitivity is $\Delta f = 1$ and the Laplace scale is simply $b = 1/\varepsilon$. If a rider can appear in many cells within one release, you must either cap contributions (keep the first $k$ points per subject) or inflate sensitivity to $k$ — otherwise the noise you add is calibrated to a lie and the guarantee is void.
• A declared release budget. Each histogram release spends $\varepsilon$ once. A dashboard that refreshes hourly composes those spends over time, so the per-release $\varepsilon$ must be drawn from a running ledger, not chosen ad hoc. Stand up the accountant described in privacy budget management and treat each grid publication as a debit against a per-subject cumulative ceiling.
• A fixed, public grid origin. The grid’s origin, cell size, and extent must be fixed in advance and independent of the data. Snapping the grid to the data’s bounding box leaks the extent of the point set; publishing a different resolution per release lets an adversary difference two releases. Pin the grid to a published reference (for example, a standard geohash or H3 resolution) so it is reproducible and auditable.

Step 1: Choose Grid Resolution from the Sensitivity Tier

Resolution is a privacy decision before it is a cartographic one. A finer grid produces more cells, each holding a smaller true count, so the fixed noise magnitude $b = 1/\varepsilon$ becomes a larger fraction of every cell — relative error scales roughly as $b / c$ for a cell of true count $c$. Halving the cell size in each dimension quarters the average count and roughly quadruples the relative error, all before you have spent a single extra unit of budget. The higher a tile’s sensitivity tier — as scored by the models feeding how to calculate spatial k-anonymity thresholds — the coarser its minimum cell must be, so that even the emptiest published cell keeps an expected population above the k-anonymity floor and the noise cannot single out a dwelling.

The mapping below is the resolution floor by tier; the pipeline may go coarser but never finer than the tier permits.

from __future__ import annotations
import math
from typing import Dict

# Minimum cell edge in decimal degrees, keyed by sensitivity tier.
# Higher tier -> coarser floor so sparse cells keep a defensible population.
# At mid-latitude, 0.01 deg latitude is ~1.1 km; 0.005 deg is ~0.55 km.
TIER_MIN_RESOLUTION_DEG: Dict[int, float] = {
    0: 0.0025,   # public context: ~275 m cells
    1: 0.005,    # low:            ~550 m cells
    2: 0.01,     # medium:         ~1.1 km cells
    3: 0.02,     # high:           ~2.2 km cells — uniform release only under a tight budget
}

def resolution_for_tier(tier: int, requested_deg: float) -> float:
    """Clamp a requested cell size up to the tier's coarseness floor.

    Returning a coarser cell than requested is always privacy-safe; returning
    a finer one would raise per-cell relative noise and risk sub-floor cells.
    """
    if tier not in TIER_MIN_RESOLUTION_DEG:
        raise ValueError(f"unknown sensitivity tier: {tier}")
    floor: float = TIER_MIN_RESOLUTION_DEG[tier]
    return max(requested_deg, floor)  # never finer than the floor

The max — not min — is the load-bearing line: a caller who asks for a 100 m grid on Tier 2 data is silently coarsened to 1.1 km, because honouring the request would publish cells whose expected count is too low to hide anyone. Record the resolution actually used alongside the release; it is audit evidence that the tier floor was respected.

Step 2: Bin Points to the Uniform Grid

With the resolution fixed, binning is a deterministic assignment of every point to the cell whose bounds contain it. numpy.histogram2d does this in one vectorised pass over the coordinate arrays, and because the grid edges are derived from the fixed origin and resolution rather than from the data’s min/max, two releases over shifting data still align cell-for-cell. The output is a dense integer count matrix — the true histogram, which never leaves the curator.

import numpy as np
import geopandas as gpd
from dataclasses import dataclass

@dataclass(frozen=True)
class GridSpec:
    """A fixed, data-independent grid. Origin and resolution are public."""
    minx: float
    miny: float
    maxx: float
    maxy: float
    resolution: float

    @property
    def x_edges(self) -> np.ndarray:
        n: int = math.ceil((self.maxx - self.minx) / self.resolution)
        return self.minx + self.resolution * np.arange(n + 1)

    @property
    def y_edges(self) -> np.ndarray:
        n: int = math.ceil((self.maxy - self.miny) / self.resolution)
        return self.miny + self.resolution * np.arange(n + 1)

def bin_points(points: gpd.GeoDataFrame, grid: GridSpec) -> np.ndarray:
    """Return the true (noise-free) count matrix for points over the grid.

    Points outside the fixed extent are dropped, not clipped to the edge,
    so an out-of-bounds coordinate cannot inflate a boundary cell.
    """
    if points.crs is None or points.crs.to_epsg() != 4326:
        raise ValueError("points must be EPSG:4326 decimal degrees")
    xs: np.ndarray = points.geometry.x.to_numpy()
    ys: np.ndarray = points.geometry.y.to_numpy()
    counts, _, _ = np.histogram2d(xs, ys, bins=[grid.x_edges, grid.y_edges])
    return counts.astype(np.int64)  # exact integer counts, curator-only

Binning is where a bounded-contribution policy is actually enforced: if a subject can contribute multiple points, cap or aggregate them before this call, because once points are anonymous marks on a grid the histogram cannot tell one busy rider from ten quiet ones. The same logic underpins geohash- or H3-cell binning — swap the rectangular edges for the cell index of your discrete global grid and the rest of the pipeline is unchanged.

Step 3: Add Per-Cell Laplace Noise Under a Per-Release ε

Every cell in the true matrix gets an independent Laplace draw at scale $b = \Delta f / \varepsilon$. Because a single subject affects at most one cell by at most one count, the whole matrix is a $\Delta f = 1$ query and one $\varepsilon$ covers the entire grid — you do not split the budget across cells. This is the parallel-composition property of histograms and it is what makes uniform-grid aggregation cheap: a million-cell grid costs the same $\varepsilon$ as a nine-cell grid. The draw itself is a standard Laplace mechanism; the calibration details and the Gaussian alternative for $\ell_2$-bounded queries live in spatial noise mechanisms.

def add_laplace_noise(
    counts: np.ndarray,
    epsilon: float,
    sensitivity: float = 1.0,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Add i.i.d. Laplace(0, sensitivity/epsilon) noise to every cell.

    One epsilon covers the whole grid: a subject touches at most one cell,
    so this is a single sensitivity-1 histogram query (parallel composition),
    NOT one query per cell. Returns floats — cells may now be negative.
    """
    if epsilon <= 0:
        raise ValueError("epsilon must be positive")
    gen: np.random.Generator = rng or np.random.default_rng()
    scale: float = sensitivity / epsilon         # b = Δf / ε
    noise: np.ndarray = gen.laplace(loc=0.0, scale=scale, size=counts.shape)
    return counts.astype(np.float64) + noise      # noised, pre-repair

The budget debit happens exactly once here, and it must be recorded before the noised matrix is handed downstream. Everything after this line is post-processing on a released quantity and is free — a fact the next step relies on completely.

Step 4: Post-Process — Clamp to ≥0 and Reconcile the Total

The noise has broken two things a consumer expects from a count surface: some cells are now negative, and the cell sum no longer matches the true total. Both repairs are pure post-processing on the DP output, so by the post-processing invariance of differential privacy they cost nothing and cannot weaken the guarantee. Clamping negatives to zero is unconditional. Total-count consistency — rescaling the non-negative cells so they sum to a separately-noised grand total — is optional and worth it when the dashboard reports both per-cell and headline figures, because an inconsistent pair (cells sum to 4,102 but total says 4,088) reads as a bug and invites the consumer to difference them.

def postprocess(
    noised: np.ndarray,
    noisy_total: float | None = None,
) -> np.ndarray:
    """Clamp to non-negative and optionally rescale to a noisy total.

    Post-processing invariance: these transforms touch only the already-DP
    matrix, so they spend zero additional budget. Clamp is always applied;
    consistency rescaling runs only when a (separately-noised) total is given.
    """
    clamped: np.ndarray = np.maximum(noised, 0.0)   # non-negativity — free
    if noisy_total is None:
        return clamped
    current: float = float(clamped.sum())
    if current <= 0.0:
        return clamped                              # nothing to rescale onto
    target: float = max(noisy_total, 0.0)
    return clamped * (target / current)             # preserve the grand total

Two cautions. First, clamping is a biased operation: it can only push the expected released count up, so a grid of mostly-empty cells develops a systematic positive bias that inflates apparent density in sparse regions — the sparse-region failure mode addressed head-on below. Second, if you enforce consistency, the total must itself be a noised quantity spent from the budget (or derived from the same noised cells), never the true total — releasing the exact grand total alongside noised cells hands the adversary a clean constraint.

Step 5: Measure Utility — Relative Error and Preserved Hotspots

A DP release is only useful if the dashboard’s decisions survive the noise, so utility is measured against the analytics task, not in the abstract. Two metrics carry most of the weight. Population-weighted relative error answers “how wrong is a typical busy cell?” and is dominated by the dense cells a planner actually acts on. Hotspot preservation — the overlap between the true top-$k$ cells and the released top-$k$ — answers “does the map still point at the right places?”, which is often the only thing a heatmap consumer cares about. Reporting both prevents the classic mistake of accepting a release whose average error looks fine while its ranking of hotspots has been scrambled by sparse-cell noise.

from typing import NamedTuple

class UtilityReport(NamedTuple):
    weighted_relative_error: float   # sum|released-true| / sum(true)
    hotspot_overlap: float           # |top-k true ∩ top-k released| / k
    max_abs_error: float

def utility_report(true_counts: np.ndarray, released: np.ndarray, k: int = 10) -> UtilityReport:
    """Compare a released grid against the (curator-only) truth.

    Population-weighted error tracks the dense cells planners act on; hotspot
    overlap tracks whether the ranking of busy cells survived the noise.
    """
    true_flat: np.ndarray = true_counts.ravel().astype(np.float64)
    rel_flat: np.ndarray = released.ravel().astype(np.float64)
    total: float = float(true_flat.sum()) or 1.0
    wre: float = float(np.abs(rel_flat - true_flat).sum() / total)
    k = min(k, true_flat.size)
    top_true = set(np.argsort(true_flat)[-k:])
    top_rel = set(np.argsort(rel_flat)[-k:])
    overlap: float = len(top_true & top_rel) / k
    return UtilityReport(wre, overlap, float(np.abs(rel_flat - true_flat).max()))

If the weighted relative error clears the analytics SLA but hotspot overlap is poor, the grid is too fine for the budget: coarsen the resolution (Step 1) so each cell holds more mass, or move to a data-adaptive partition. That second option — spending a slice of the budget to learn where resolution should be fine and where it should be coarse — is exactly the private quadtree decomposition for spatial histograms technique, and it is the standard escape hatch when a single uniform resolution cannot serve both downtown and the exurbs.

Step 6: Assemble and Validate the Release

The pieces compose into one PrivateSpatialHistogram that takes raw points and a budget and returns a GeoDataFrame of noised cell polygons ready for the dashboard, with an internal _validate() that fails closed on any broken invariant. The orchestration is deliberately linear — bin, noise, repair, package — so the single budget debit is easy to audit and impossible to double-spend.

class PrivateSpatialHistogram:
    """Produce an ε-DP density grid from raw points via a uniform partition."""

    def __init__(self, grid: GridSpec, epsilon: float, sensitivity: float = 1.0):
        if epsilon <= 0 or sensitivity <= 0:
            raise ValueError("epsilon and sensitivity must be positive")
        self.grid = grid
        self.epsilon = epsilon
        self.sensitivity = sensitivity
        self._spent: float = 0.0

    def release(
        self,
        points: gpd.GeoDataFrame,
        enforce_consistency: bool = True,
        rng: np.random.Generator | None = None,
    ) -> gpd.GeoDataFrame:
        """Bin, noise once, repair, and emit noised cell polygons (EPSG:4326)."""
        gen = rng or np.random.default_rng()
        true_counts: np.ndarray = bin_points(points, self.grid)
        noised: np.ndarray = add_laplace_noise(
            true_counts, self.epsilon, self.sensitivity, gen)
        self._spent += self.epsilon                 # single debit for the grid
        noisy_total: float | None = None
        if enforce_consistency:                     # spend nothing extra: reuse cells
            noisy_total = float(np.maximum(noised, 0.0).sum())
        released: np.ndarray = postprocess(noised, noisy_total)
        return self._to_geoframe(released, points.crs)

    def _to_geoframe(self, matrix: np.ndarray, crs) -> gpd.GeoDataFrame:
        from shapely.geometry import box
        res = self.grid.resolution
        recs = []
        for ix in range(matrix.shape[0]):
            for iy in range(matrix.shape[1]):
                x0 = self.grid.minx + ix * res
                y0 = self.grid.miny + iy * res
                recs.append({"count": float(matrix[ix, iy]),
                             "geometry": box(x0, y0, x0 + res, y0 + res)})
        return gpd.GeoDataFrame(recs, geometry="geometry", crs=crs)

    def _validate(self) -> None:
        """Fail closed on any broken DP or grid invariant. Safe in CI."""
        gen = np.random.default_rng(7)
        pts = gpd.GeoDataFrame(
            {"geometry": gpd.points_from_xy(
                gen.uniform(-0.02, 0.02, 4000), gen.uniform(-0.02, 0.02, 4000))},
            crs="EPSG:4326")
        out = self.release(pts, enforce_consistency=True, rng=gen)
        assert (out["count"] >= 0).all(), "non-negativity violated post-clamp"
        assert math.isclose(self._spent, self.epsilon, abs_tol=1e-9), "budget mis-debited"
        assert resolution_for_tier(2, 0.001) == 0.01, "tier floor not enforced"
        # noise is unbiased pre-clamp: mean released total tracks true total loosely
        assert out["count"].sum() > 0, "release collapsed to empty grid"
        print("PrivateSpatialHistogram invariants hold")


if __name__ == "__main__":
    spec = GridSpec(minx=-0.02, miny=-0.02, maxx=0.02, maxy=0.02, resolution=0.01)
    PrivateSpatialHistogram(spec, epsilon=1.0)._validate()

Threat Model Considerations

Uniform-grid aggregation defends against a curator-trusted adversary who sees only the released grid and its metadata. The surfaces specific to spatial histograms:

• Difference-of-releases attacks. Two hourly grids over the same population, each $\varepsilon$-DP, compose to $2\varepsilon$; an adversary who differences them recovers movement. The defence is composition accounting against a per-subject ceiling via privacy budget management, not per-release reasoning.
• Data-dependent grid leakage. Snapping the grid origin or extent to the data reveals the bounding box of the point set — a free, un-noised signal. A fixed public grid closes this; it is why Step 1 forbids finer-than-floor and Step 2 derives edges from a fixed origin.
• Sparse-cell singling-out. In a near-empty cell, a released count of 1 against an expected 0 can betray a lone subject even after noise, because the clamp bias makes small positive values suspicious. The k-anonymity floor from how to calculate spatial k-anonymity thresholds sets the minimum cell size that keeps expected populations above the singling-out threshold.
• Consistency-constraint exploitation. Publishing the true grand total alongside noised cells gives the adversary an exact linear constraint that reduces the effective noise. Any total used for reconciliation must itself be noised or derived from the noised cells.
• Unbounded contribution. If a subject can land in many cells and you calibrated to $\Delta f = 1$, the real sensitivity is higher and the guarantee is broken silently. Contribution capping at binning time is the only fix.

Validation & Compliance Checklist

Wire each control into CI with a measurable pass criterion rather than reviewing by eye:

1. Non-negativity — PASS if every released cell is $\ge 0$ across a $\ge 10^4$-cell property test after postprocess. Clamping is unconditional; fail the build on any negative.
2. Single-debit budgeting — PASS if one release() call increments spend by exactly $\varepsilon$ regardless of cell count, proving parallel composition is used rather than per-cell splitting.
3. Tier-floor enforcement — PASS if resolution_for_tier(tier, r) never returns a value below the tier’s floor for any requested r, so no release publishes sub-floor cells.
4. Data-independent grid — PASS if the grid edges are identical across two releases over different point sets with the same GridSpec, confirming the partition leaks nothing about the data extent.
5. Consistency safety — PASS if the total used for reconciliation is a noised or cell-derived quantity, never the true sum; assert the reconciliation input is not true_counts.sum().
6. Utility SLA — PASS if population-weighted relative error and top-$k$ hotspot overlap both clear the documented analytics thresholds; reject the parameter set otherwise.

For regulated releases, bind these to a clause: GDPR Article 25 data-minimisation caps the per-subject cumulative $\varepsilon$, and any HIPAA-adjacent health-mobility release must keep the smallest published geography above the Safe Harbor 20,000-population floor — a constraint that resolves directly to a minimum cell size in Step 1.

Failure Modes & Remediation

Aggregate DP fails quietly — it returns a plausible-looking surface that is either useless or unsafe. The high-frequency production failures:

• Sparse-region noise dominance. Suburban and rural cells hold true counts near zero, so noise at scale $1/\varepsilon$ plus clamp bias manufactures phantom hotspots in empty land. Detection: hotspot overlap collapses while dense-cell error looks fine. Remediation: coarsen the resolution for low-density extents, or switch to the adaptive private quadtree decomposition so budget concentrates where data lives.
• Over-fine grid burning relative accuracy. A cartographically appealing fine grid quarters cell counts and multiplies relative error without spending more $\varepsilon$. Detection: weighted relative error rises sharply as resolution drops below the tier floor. Remediation: honour the Step 1 floor; treat resolution as a privacy parameter, not a rendering preference.
• Silent sensitivity underestimation. A subject contributes to many cells but noise was calibrated to $\Delta f = 1$. Detection: audit the maximum cells-per-subject at binning time. Remediation: cap contributions or raise sensitivity to the observed maximum and re-noise.
• Budget exhaustion under dashboard refresh. An hourly heatmap composes $\varepsilon$ every hour and silently blows the per-subject ceiling within a day. Detection: the accountant crosses its cap. Remediation: lengthen the refresh interval, coarsen releases to spend less, or rotate to a fresh cohort — see privacy budget management.
• Inconsistent headline and cells. The dashboard’s grand total and its cell sum disagree, inviting a differencing bug report or attack. Detection: released total minus cell sum is non-zero. Remediation: enable consistency reconciliation against a noised total.

Frequently Asked Questions

Related

This guide is part of the Differential Privacy for Geospatial Data section — start there for how noise mechanisms, aggregation, geo-indistinguishability, and trajectory privacy fit together.

• Building Differentially Private Heatmaps — turning the noised grid produced here into a rendered, interpolated surface.
• Private Quadtree Decomposition for Spatial Histograms — the data-adaptive alternative when one uniform resolution cannot serve dense and sparse regions at once.
• Spatial Noise Mechanisms — how the Laplace and Gaussian noise added per cell is calibrated.
• Privacy Budget Management — composing epsilon across repeated dashboard releases against a per-subject ceiling.
• How to Calculate Spatial K-Anonymity Thresholds — the cell-population floor that sets the coarsest safe resolution.

Up: Differential Privacy for Geospatial Data