Spatial Epi Public Health GIS

Adaptive Areal Aggregation for Rare-Disease Counts

Fixed census-tract reporting leaves most rare-disease cells below a publishable threshold, so a region-growing merge that absorbs contiguous neighbors until each region meets both a minimum population and a minimum case count is the defensible alternative. This guide, part of spatial k-anonymity for health microdata, builds that aggregation with libpysal contiguity weights, keeps every region connected, and orders the whole procedure deterministically so a reviewer can reproduce it exactly.

Problem Context & Constraints

Report a rare disease at fixed census-tract resolution and the map is mostly holes. A condition with a national incidence of a few per hundred thousand produces tracts of zero, one, or two cases almost everywhere; publish those counts and you have released a near-unique record for every non-zero tract. Suppressing all the small cells — the reflex response — turns the map into Swiss cheese and destroys the very spatial signal the release was meant to carry. Neither raw publication nor blanket suppression is acceptable.

The joint nature of the constraint is what makes rare-disease aggregation distinct from ordinary small-area work:

  • Population alone is not enough. A region can hold fifty thousand residents and still contain exactly one case. A lone case in a large region is still a potential re-identification when the condition is distinctive, so the case count is an independent gate that population size cannot satisfy.
  • Case count alone is not enough. A region can accumulate five cases while covering so few residents that each case is nearly identifiable. Both floors — minimum at-risk population for k-anonymity and minimum case count for cell safety — must hold simultaneously.
  • Contiguity is mandatory. Merging a sparse tract with a distant populous one to hit the thresholds fabricates a geography no epidemiologist can interpret; regions must be spatially connected to remain meaningful on a map and defensible in review.

The tool for this is region-growing — the same family as the Max-P regionalization problem, which partitions areal units into the largest number of contiguous regions each satisfying a floor constraint. A full Max-P solve is combinatorial; in production a greedy contiguity merge driven by spatial weights matrix construction gives a deterministic, auditable approximation that clears both floors while keeping regions connected. It pairs naturally with dasymetric population suppression: suppression removes cells that cannot be published, and adaptive aggregation grows what remains into publishable regions.

Region Growing to a Joint Threshold Three stages left to right. Stage one is a single tract with one case and 380 residents, below threshold. Stage two adds one contiguous neighbor, reaching three cases and 1,200 residents, still below threshold. Stage three has grown to a two-by-two block of four contiguous tracts totaling five cases and 5,400 residents, which meets both the minimum case count of five and the minimum population, so growth stops. A region grows until it clears both floors 1 cases 1 · pop 380 below threshold +neighbor 1 2 cases 3 · pop 1,200 still below +neighbors 1 2 1 1 cases 5 · pop 5,400 meets min_cases & min_pop Growth absorbs contiguous neighbors deterministically, then stops the moment both floors are met

Prerequisites

Pin the stack so region labels and merge order are reproducible across reviewers:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x, pyproj 3.6.x)
  • libpysal 4.12.1 — contiguity weights
  • numpy 1.26.4

Input state assumed below:

  • Tract polygons with a stable identifier (geoid), an at-risk population field, and a rare-disease case count, all in one projected CRS so any area or adjacency logic is geometrically sound.
  • Valid, single-typed geometries with clean shared boundaries; a sliver that splits an edge drops a real adjacency and can strand a tract that should have merged.
  • Two written thresholds — min_pop for k-anonymity and min_cases for cell safety — plus a rule that ties break on the stable identifier so the result is deterministic.

Step-by-Step Solution

Build Queen contiguity from the tracts, then grow regions greedily: process seed tracts in ascending identifier order, and while a seed’s region falls short of either floor, absorb the contiguous neighboring region that best strengthens it, breaking ties on the identifier. Because both the seed order and the tie-break are fixed, the same input always yields the same regionalization.

# Deterministic region-growing merge to a joint population + case-count floor.
# Pinned: geopandas==0.14.4, libpysal==4.12.1, numpy==1.26.4
import logging
import geopandas as gpd
from libpysal.weights import Queen

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

def grow_regions(tracts, *, min_pop, min_cases, geoid="geoid",
                 pop_field="pop_at_risk", case_field="cases"):
    """
    Grow contiguous regions until each meets BOTH min_pop and min_cases.
    Deterministic: seeds processed by ascending geoid; ties break on geoid.
    Returns (labels_by_geoid, merge_log).
    """
    tracts = tracts.sort_values(geoid).reset_index(drop=True)
    ids = list(tracts[geoid].astype(str))
    w = Queen.from_dataframe(tracts, ids=ids)              # contiguity graph

    pop = dict(zip(ids, tracts[pop_field].astype(float)))
    cases = dict(zip(ids, tracts[case_field].astype(float)))
    region = {u: u for u in ids}                            # each tract seeds its own region
    members = {u: [u] for u in ids}
    merge_log = []

    def stats(r):
        return (sum(pop[m] for m in members[r]),
                sum(cases[m] for m in members[r]))

    def meets(r):
        p, c = stats(r)
        return p >= min_pop and c >= min_cases

    for seed in ids:                                        # ascending, stable order
        r = region[seed]
        while not meets(r):
            # regions contiguous to r (across any member boundary), excluding r itself
            adj = set()
            for m in members[r]:
                for nb in w.neighbors[m]:
                    if region[nb] != r:
                        adj.add(region[nb])
            if not adj:
                log.warning("region %s isolated below threshold; flag for review", r)
                merge_log.append({"region": r, "action": "isolated_below_threshold"})
                break
            # prefer the weakest contiguous region (lowest pop), tie-break on id
            target = min(adj, key=lambda rr: (stats(rr)[0], rr))
            # fold target into r, keeping the lexicographically smaller label as canonical
            keep, drop = sorted((r, target))
            members[keep].extend(members[drop])
            for m in members[drop]:
                region[m] = keep
            members.pop(drop, None)
            r = keep
            merge_log.append({"kept": keep, "absorbed": drop,
                              "pop": stats(keep)[0], "cases": stats(keep)[1]})
            log.info("region %s absorbed %s -> pop=%.0f cases=%.0f",
                     keep, drop, *stats(keep))

    log.info("regionalization complete: regions=%d from tracts=%d",
             len(set(region.values())), len(ids))
    return region, merge_log

def dissolve_to_regions(tracts, region, geoid="geoid",
                        pop_field="pop_at_risk", case_field="cases"):
    """Materialize region polygons with summed population and cases."""
    tracts = tracts.copy()
    tracts["region"] = tracts[geoid].astype(str).map(region)
    regions = tracts.dissolve(by="region",
                              aggfunc={pop_field: "sum", case_field: "sum"})
    return regions.reset_index()

The greedy rule folds each undersized region into its weakest contiguous neighbor, which consolidates two fragile areas rather than diluting a healthy one, and it always keeps the lexicographically smaller identifier as the region label so the output naming is itself deterministic. A region with no contiguous neighbor — an island, or a tract fully enclosed by already-suppressed cells — is flagged rather than force-merged, because inventing a non-adjacent link would produce a region that does not exist on the ground.

Validation & Edge Cases

Three assertions must pass before any regionalized layer is released.

1. Every region meets both floors. This is the whole point of the exercise, so it is a release-blocking check, not a warning:

def assert_thresholds(regions, min_pop, min_cases,
                      pop_field="pop_at_risk", case_field="cases"):
    bad = regions[(regions[pop_field] < min_pop) | (regions[case_field] < min_cases)]
    assert bad.empty, f"{len(bad)} regions below a floor; do not release"
    return True

If bad is non-empty it is almost always an isolated region flagged during growth; those must be resolved by policy — merged across a non-contiguity exception with explicit sign-off, or suppressed entirely — never quietly published.

2. Contiguity is preserved. Each dissolved region must be a single connected polygon. A MultiPolygon with disjoint parts means a merge jumped a gap and the region is not actually connected:

from shapely.geometry import MultiPolygon

def assert_contiguous(regions):
    for _, row in regions.iterrows():
        g = row.geometry
        parts = len(g.geoms) if isinstance(g, MultiPolygon) else 1
        assert parts == 1, f"region {row['region']} is disconnected ({parts} parts)"
    return True
INFO adaptive_agg: regionalization complete: regions=142 from tracts=1174

A legitimate exception is genuine archipelago geography, where a single administrative region really does span islands; handle those with an explicit allow-list keyed on the identifier rather than by loosening the global check.

3. Deterministic ordering holds. Re-run the growth on the same input and the region mapping must be identical. Because seeds are processed in ascending geoid order and ties break on the identifier, two runs — or two analysts — produce byte-identical labels. Sorting the input on geoid before anything else is what guarantees this; skip it and the Queen graph’s neighbor order can drift and quietly change which weak neighbor gets absorbed first.

A subtler failure mode is threshold interaction: setting min_cases high relative to disease rarity forces enormous regions that swallow whole counties, erasing spatial signal. If the regionalization collapses to a handful of giant regions, the thresholds are fighting the data — revisit whether the case-count floor is proportionate to the condition’s incidence before accepting a map that says almost nothing.

Compliance Notes

Three artifacts make an adaptive aggregation defensible under review:

  • The thresholds, logged with the release. Record min_pop and min_cases in the audit block; a region map is uninterpretable — and potentially disclosive — if a later consumer does not know which floors it was built to clear. Log them as structured fields, not prose.
  • The ordered merge log. Because the greedy growth is order-dependent, the sequence of absorptions is the method. Persisting merge_log alongside the pinned seed-ordering rule lets a reviewer replay the exact regionalization; without it the output cannot be reproduced.
  • Isolated-region disposition. Every region flagged as isolated-below-threshold must have a recorded resolution — merged by exception, or suppressed — so the file carries no silently under-protected region. That disposition list is the auditable record of the hard cases the automation could not resolve on its own.

Bind the thresholds and library versions to the output with a configuration hash, and treat the threshold-passing assertion as a gate that halts the pipeline on failure rather than a check that merely logs. A rare-disease map is one of the highest-risk releases a program produces; the audit trail is what lets it withstand the scrutiny that risk invites.