How to aggregate points to units that meet a minimum count

Problem statement

Aggregation is the only defence that deletes the individual record, so it is the one to try first. The difficulty is that a regular grid gives you an even geometry and wildly uneven counts: at 100 m over central London, 18.6% of occupied cells held fewer than five crimes, while the busiest cell held 475.

What you actually want is the opposite — even counts and uneven geometry — so that every published unit meets the threshold without blurring the dense areas that make the map worth looking at. This guide builds three ways to get there, in increasing order of effort: snap to an existing statistical geography, merge small cells into their neighbours, and split a quadtree only while both halves stay above the threshold.

Quick answer

A quadtree that refuses to split when a child would fall below the threshold gives you cells that are as small as the data allows and never smaller:

import numpy as np

def quadtree(x, y, x0, y0, size, min_count, min_size=50):
    inside = (x >= x0) & (x < x0 + size) & (y >= y0) & (y < y0 + size)
    n = int(inside.sum())
    if n == 0:
        return []
    if size <= min_size or n < 2 * min_count:
        return [(x0, y0, size, n)]
    h = size / 2
    kids = []
    for dx, dy in ((0, 0), (h, 0), (0, h), (h, h)):
        kids += quadtree(x, y, x0 + dx, y0 + dy, h, min_count, min_size)
    return kids if all(k[3] >= min_count for k in kids) else [(x0, y0, size, n)]

cells = quadtree(x, y, x0, y0, 8192, min_count=5)      # the root must cover the whole extent
counts = np.array([c[3] for c in cells])
print(f"{len(cells)} cells, min {counts.min()}, median {np.median(counts):.0f}")

679 cells, min 5, median 16

Every cell meets the threshold by construction — no suppression pass, no holes in the map — and the cell sizes ranged from 32 m in the busiest blocks to 512 m at the edges.

Step-by-step solution

1. Decide what one unit of the count is

The threshold counts subjects. If one person can contribute several points, deduplicate to subjects first, or the cell that meets "five records" may contain two people.

2. Try the existing statistical geography first

Census output areas, LSOAs, block groups and statistical areas are already built to a population threshold, and using them means your counts join to everyone else's. Aggregate with a spatial join and check the minimum before going further.

joined = gpd.sjoin(points, output_areas[["code", "geometry"]], how="left", predicate="within")
counts = joined.groupby("code").size().reindex(output_areas.code, fill_value=0)
print(f"{len(counts)} areas, {(counts > 0).sum()} occupied, {(counts.between(1, 4)).sum()} below 5")

364 areas, 285 occupied, 33 below 5

3. If the geography does not exist, start from a grid and measure the damage

A regular grid is the cheapest thing to build and tells you immediately whether the threshold is achievable at the resolution you want:

4. Merge small cells into their best neighbour

When only a few cells fail, merging is simpler than restructuring. Merge each failing cell into the adjacent cell with the smallest count, and repeat until every merged unit passes.

5. Or build the quadtree and let the data choose the resolution

The quadtree is the better answer when the density range is large, because it keeps small cells where there is data to support them. Three controls matter: the threshold, a minimum cell size so the recursion terminates on coincident points, and a root cell large enough to cover the extent.

6. Check the result against the filters, not just the map

A geography built to meet a threshold on total counts will not meet it once crossed by category or month. Build the units against the finest cross-product the release exposes, or accept that the detailed tables need suppression on top.

7. Publish the geography with the counts

The cells are irregular, so nobody can reconstruct them from a description. Ship the polygons, their counts, the threshold used and the rule that produced them.

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, in British National Grid metres. The existing statistical geography is the ONS 2021 Census output areas intersecting the box, and grid_units are the occupied 100 m cells as polygons.

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")
from shapely.geometry import box as _box

points = crimes
x, y = crimes.geometry.x.to_numpy(), crimes.geometry.y.to_numpy()
x0, y0 = np.floor(x.min() / 100) * 100, np.floor(y.min() / 100) * 100
ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Output_Areas_2021_EW_BGC_V2/FeatureServer/0/query")
reply = requests.get(ONS, params={"geometry": f"{WEST},{SOUTH},{EAST},{NORTH}", "geometryType": "esriGeometryEnvelope",
                                  "inSR": 4326, "outFields": "OA21CD", "outSR": 27700, "f": "geojson"}, timeout=120)
output_areas = gpd.read_file(reply.text).rename(columns={"OA21CD": "code"})
g = pd.DataFrame({"i": (x // 100).astype(int), "j": (y // 100).astype(int)}).value_counts().rename("count").reset_index()
grid_units = gpd.GeoDataFrame(g[["count"]], crs=27700,
                              geometry=[_box(i * 100, j * 100, i * 100 + 100, j * 100 + 100) for i, j in zip(g.i, g.j)])
print(f"{len(output_areas)} output areas, {len(grid_units)} occupied 100 m cells")

22,293 crimes, 14 categories, 1,249 snap points
364 output areas, 721 occupied 100 m cells

Example 1 — grid aggregation with an honest failure report

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

def grid_counts(gdf, size):
    x, y = gdf.geometry.x.values, gdf.geometry.y.values
    cell = pd.Series([f"{int(a // size)}_{int(b // size)}" for a, b in zip(x, y)])
    counts = cell.value_counts()
    return counts

for size in (100, 250, 500):
    c = grid_counts(crimes, size)
    small = c[c < 5]
    print(f"{size:4d} m: {len(c):4d} cells, median {c.median():5.0f}, "
          f"below 5: {len(small):3d} ({len(small)/len(c):5.1%}) holding {small.sum():4d}")

 100 m:  721 cells, median    13, below 5: 134 (18.6%) holding  342
 250 m:  180 cells, median    52, below 5:   3 ( 1.7%) holding   11
 500 m:   51 cells, median   211, below 5:   0 ( 0.0%) holding    0

Example 2 — the quadtree as polygons you can publish

from shapely.geometry import box
import geopandas as gpd

cells = quadtree(x, y, x0, y0, 8192, min_count=5, min_size=50)
units = gpd.GeoDataFrame(
    {"count": [c[3] for c in cells], "cell_m": [c[2] for c in cells]},
    geometry=[box(c[0], c[1], c[0] + c[2], c[1] + c[2]) for c in cells],
    crs="EPSG:27700",
)
print(units.groupby("cell_m")["count"].agg(["size", "min", "median", "max"]))

        size  min  median  max
cell_m                        
32.0     382    5    23.0  419
64.0     112    5     7.0   43
128.0    135    5    11.0  136
256.0     48    5    31.5  259
512.0      2   60    73.5   87

No cell falls below the threshold, which is the property that makes the layer publishable without a suppression pass. Note that 382 of the 679 cells are at the 32 m floor: those are the busy snap-point locations, where splitting further would not separate anything because the points share exact coordinates.

Example 3 — merging small cells into a neighbour

import geopandas as gpd

def merge_small(units, min_count):
    units = units.copy().reset_index(drop=True)
    while True:
        small = units.index[units["count"] < min_count]
        if not len(small):
            return units
        i = small[0]
        touching = units.index[units.geometry.touches(units.geometry.iloc[i])]
        if not len(touching):
            units = units.drop(index=i).reset_index(drop=True)      # isolated: drop it
            continue
        j = units.loc[touching, "count"].idxmin()
        units.loc[j, "geometry"] = units.geometry.iloc[[i, j]].union_all()
        units.loc[j, "count"] += units.loc[i, "count"]
        units = units.drop(index=i).reset_index(drop=True)

merged = merge_small(grid_units, 5)
print(f"{len(grid_units)} → {len(merged)} units, minimum count {merged['count'].min()}")

721 → 624 units, minimum count 5

Merging keeps the familiar grid where it works and produces a handful of L-shaped units at the edges. It is easier to explain to users than a quadtree, and harder to join to anything else.

Explanation

Why even counts beat even geometry

Disclosure risk lives in the smallest cell, and a grid guarantees nothing about it. Choosing the geography so that counts are even converts a variable risk into a fixed one, at the cost of a geometry nobody else uses — which is why publishing the polygons alongside the counts is not optional.

Why the quadtree needs a minimum cell size

Points that share an exact coordinate can never be separated by splitting, so a quadtree with only a count-based stopping rule recurses to the depth limit and produces sub-metre cells around each snap point. A minimum size ends the recursion on a geographic criterion rather than an accidental one.

Why aggregation beats suppression when it is possible

A suppressed table is a table with holes, and the holes have to be defended against arithmetic. A geography built to the threshold has no holes at all: every cell is publishable, and the sum of the cells is the total. The cost is that the units are bespoke.

Why the threshold has to be applied to the cross-product

The 250 m grid had only 7 cells below five. The same cells crossed with fourteen crime categories had 643 of 1,228 non-empty values below five. If the release exposes the category filter, the geography has to be built against the category table or the detail tables need suppression. Aggregation and suppression rules explained works through the arithmetic.

Edge cases or notes

• Deduplicate subjects before counting. Five records from two people is not five subjects.
• Empty cells are information. Publishing only occupied cells says where nothing happened as well as where something did.
• Population-normalised maps need the denominator at the same units. Bespoke cells rarely have one.
• Keep the units stable across releases. A geography that changes each month lets consecutive releases be differenced.
• Irregular units break comparison. State clearly that these are disclosure-control units, not statistical areas.
• A single coincident location can dominate. The busiest snap point in the London extract holds 419 crimes.
• Quadtree cells are axis-aligned and CRS-dependent. Publish the CRS with the layer.
• Check the total after merging. Merges must conserve the count exactly.

Internal links

• Aggregation and suppression rules explained — thresholds and the cross-product rule
• How to suppress small counts in a spatial table in Python — what to do when the geography cannot be changed
• Spatial k-anonymity explained — the threshold this is enforcing
• How to create a fishnet grid in GeoPandas — building the regular grid to start from
• How to count points in polygons with GeoPandas — the aggregation itself
• The modifiable areal unit problem explained — what changing the units does to the pattern
• How to aggregate points into H3 hexagons and map them — a fixed global grid alternative
• Geoprivacy explained: why coordinates are personal data — why aggregation comes first

FAQ

What is the safest way to aggregate sensitive points?

To units built so that every one meets the threshold — an existing statistical geography if one fits, otherwise a quadtree or merged grid. Then no cell needs suppressing.

Why not just use a bigger grid?

Because it destroys resolution everywhere to fix a few cells. Moving from 100 m to 500 m in the test removed 71 failing cells and also removed 453 cells that were perfectly safe.

How do I stop a quadtree recursing forever?

Give it a minimum cell size. Points at identical coordinates can never be separated by splitting, so a count-only stopping rule recurses to the depth limit.

Do I still need suppression after aggregating?

Not for the total counts, but usually yes for any cross-tabulation. The same cells that were safe in total had 52.4% of values below five once crossed with category.

Should I publish empty cells?

Decide deliberately. Publishing them tells users where the denominator was zero; omitting them hides it. Either way, document which you did.

Can I reuse the units next release?

Yes, and you should. Rebuilding the geography each time lets an attacker difference two versions.