Spatial Epi Public Health GIS

Areal Interpolation & Boundary Harmonization

Health data arrives on the geography whoever collected it happened to use, and almost never on the geography the analysis needs. This guide is part of Spatial Epidemiology Fundamentals & Data Standards, and it covers how to move counts and rates between incompatible areal units, how to reconcile the same nominal geography across vintages, and how to state the error each transfer introduces instead of losing it.

Concept & Epidemiological Alignment

Areal interpolation is the transfer of an attribute from a set of source zones, where it was measured, to a set of target zones, where it is needed. Cases reported by ZIP code and denominators published by census tract is the canonical example; a 2010-vintage boundary file joined to 2020 population is another, and it is more dangerous because the join succeeds.

Every method answers the same question — how much of the source zone’s total belongs to each piece of it — and they differ only in what evidence they use to answer it. Simple areal weighting uses land area, which assumes population is uniform inside the source zone. Population-weighted interpolation uses a finer population layer, usually blocks, which assumes the finer layer is right. Dasymetric methods use ancillary land-use or building data to redistribute before weighting, which assumes the ancillary layer is current and correctly classified. Accuracy rises across that list, and so does the number of assumptions that can be wrong.

Three constraints govern whether any of them is epidemiologically defensible:

  • Counts are interpolable; rates are not. Transfer the numerator and the denominator separately and divide afterwards. Interpolating a rate directly weights each source zone equally regardless of its population, which is wrong in the specific direction that makes small rural zones dominate.
  • Interpolation cannot create information. A count moved from a ZIP to five tracts carries no more spatial detail than it had; it has been spread according to an assumption. The tract-level numbers that result look precise and are not, and any downstream cluster test will treat them as though they were measured.
  • The transfer must conserve the total. The sum of the target values has to equal the source total exactly, or the interpolation has invented or destroyed cases. This is the same conservation discipline that governs Dasymetric Population Suppression, and it is checkable in one line.
Area Weighting Against Population Weighting on One Zone A source zone holding 80 cases overlaps four target tracts. Under area weighting, the cases are split in proportion to overlapping land area, giving 30, 22, 18 and 10. Under population weighting using block-level population, the same overlap gives 8, 12, 44 and 16, because most of the population lives in the third tract while most of the land lies in the first. Both distributions sum to 80; only one of them is about people. Same 80 cases, same overlap, two very different answers by land area by block population 30 38% of area 22 27% 18 23% 10 12% 8 10% of people 12 15% 44 55% 16 20% Σ = 80 Σ = 80 The top-left tract is mostly farmland; area weighting hands it nearly four times its share of cases Conservation holds in both panels, which is why conservation alone is not a validity check

That last line is worth dwelling on. Conservation is necessary and nowhere near sufficient: an interpolation can distribute a total perfectly and still be entirely wrong about where the people are. The conservation check catches arithmetic errors; only the choice of weighting layer addresses the substantive question.

Method-Selection Table

Method Weight Assumes Use when
Areal weighting overlap area population is uniform in the source zone source zones are small and homogeneous
Population weighting block or block-group population the finer population layer is accurate a census population layer is available — the usual default
Dasymetric population restricted to habitable land land-use classification is current large source zones with substantial uninhabited area
Target-density weighting a correlated variable on the target zones the correlate genuinely predicts the attribute housing units or road density are available and the attribute tracks them
Point reallocation the underlying point locations you still hold the points the aggregation was yours to begin with — always prefer this

The last row is the one most often forgotten. If the agency holds the geocoded points, no interpolation is needed at all: re-aggregate the points to the target geography and the transfer is exact. Interpolation is a repair for data received already aggregated, and reaching for it when the points are on disk imports error for no reason.

Spatial Data Prerequisites

  • Both geographies in one metric CRS, since every method computes areas. Follow Coordinate Reference Systems for Public Health and use an equal-area projection so the weights are not latitude-dependent.
  • Valid, non-overlapping source zones. Overlapping sources double count; slivers between them lose cases. Run make_valid and check that the source union covers the target union.
  • A weighting layer whose vintage matches the source data. A 2020 block population layer applied to 2015 counts assigns cases to subdivisions that did not exist yet.
  • A stable identifier on both sides, so the crosswalk is reproducible and joinable rather than positional.

Production Implementation

# Population-weighted areal interpolation with a conservation gate.
# Pinned: geopandas==1.0.1, pandas==2.2.2, numpy==1.26.4, shapely==2.0.6
import logging
import numpy as np
import pandas as pd
import geopandas as gpd

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

EQUAL_AREA = "EPSG:5070"


def build_crosswalk(source: gpd.GeoDataFrame, target: gpd.GeoDataFrame,
                    weights: gpd.GeoDataFrame, weight_col: str,
                    src_id: str = "src_id", tgt_id: str = "tgt_id") -> pd.DataFrame:
    """Fraction of each source zone's weight that falls in each target zone.

    `weights` is a fine population layer (blocks). Its centroids are assigned to
    source and target zones, and the crosswalk is the cross-tabulation of weight.
    Using centroids rather than areal overlap keeps each block whole, which avoids
    splitting an indivisible population unit across a boundary."""
    for name, gdf in (("source", source), ("target", target), ("weights", weights)):
        if gdf.crs is None:
            raise ValueError(f"{name} has no CRS; refusing to guess")
    s = source.to_crs(EQUAL_AREA)
    t = target.to_crs(EQUAL_AREA)
    w = weights.to_crs(EQUAL_AREA).copy()
    w["geometry"] = w.geometry.representative_point()

    w = gpd.sjoin(w, s[[src_id, "geometry"]], how="left", predicate="within").drop(columns="index_right")
    w = gpd.sjoin(w, t[[tgt_id, "geometry"]], how="left", predicate="within").drop(columns="index_right")

    orphan = w[src_id].isna().sum()
    if orphan:
        log.warning("%d weight units fell outside every source zone", int(orphan))

    cw = (w.dropna(subset=[src_id, tgt_id])
            .groupby([src_id, tgt_id])[weight_col].sum().reset_index())
    totals = cw.groupby(src_id)[weight_col].transform("sum")
    cw["frac"] = np.where(totals > 0, cw[weight_col] / totals, 0.0)
    log.info("crosswalk: %d source zones -> %d target zones, %d pairs",
             cw[src_id].nunique(), cw[tgt_id].nunique(), len(cw))
    return cw[[src_id, tgt_id, "frac"]]


def interpolate_counts(counts: pd.DataFrame, crosswalk: pd.DataFrame,
                       value_cols: list[str], src_id: str = "src_id",
                       tgt_id: str = "tgt_id", tol: float = 1e-6) -> pd.DataFrame:
    """Apply the crosswalk to one or more COUNT columns, then verify conservation."""
    m = crosswalk.merge(counts, on=src_id, how="left")
    for c in value_cols:
        m[c] = m[c].fillna(0.0) * m["frac"]
    out = m.groupby(tgt_id)[value_cols].sum().reset_index()

    for c in value_cols:
        src_total, tgt_total = counts[c].sum(), out[c].sum()
        drift = abs(src_total - tgt_total)
        log.info("%s: source %.4f -> target %.4f (drift %.2e)", c, src_total, tgt_total, drift)
        if drift > tol * max(1.0, abs(src_total)):
            raise AssertionError(
                f"conservation failed for {c}: {src_total} != {tgt_total}. "
                "A source zone has no weight, or the source union does not cover the target.")
    return out

The raise on conservation failure is deliberate. The most common cause is a source zone containing no population blocks — an industrial ZIP, a water-only polygon — whose cases therefore have nowhere to go. Silently returning a smaller total turns that into a rate error nobody sees; failing loudly turns it into a five-minute investigation.

Parameter Selection & Tuning

  • Choose the weighting layer by what the attribute is about. Cases of a childhood condition are better weighted by under-18 population than by total population; occupational exposures by workplace counts rather than residence. The default of total population is a default, not a rule.
  • Prefer centroid assignment to areal overlap for the weighting layer. Splitting a census block across a boundary distributes an indivisible population unit, and block centroids are published precisely so that this is unnecessary.
  • Set the conservation tolerance in absolute terms for counts and relative terms for continuous quantities. Counts should conserve to floating-point noise; anything larger is a bug.
  • Decide whether to round. Fractional cases are correct for further arithmetic and absurd in a published table. Round at the last step, using a largest-remainder method so the rounded values still sum to the original total.

Edge Cases & Failure Modes

Source zones with zero weight. A ZIP consisting of an airport has no residential population and may still report cases — staff, travellers. Population weighting cannot place them. Either fall back to area weighting for those zones, with a flag, or exclude them and report the exclusion.

Target zones straddling the source union boundary. A target tract half outside the study area receives only the cases from the covered half, and its rate is computed against its whole population. Clip the targets to the source union, or compute rates only on fully covered targets.

Nested versus crossing geographies. ZIP codes and census tracts cross; counties and tracts nest. Where the relationship is nesting, interpolation reduces to aggregation and is exact, so check for nesting before reaching for a crosswalk.

Small numbers. Interpolating a count of 3 across five targets produces five fractional values between 0.2 and 1.4, none of which is a case. The result is arithmetically fine and epidemiologically meaningless, and it will interact badly with the suppression rules in Small-Count Cell Suppression in Rate Maps.

Error Grows With the Size Ratio, Not the Method Median absolute error in interpolated tract counts plotted against the ratio of source zone area to target zone area, for three methods. All three are near zero when source and target are similar in size. Area weighting degrades fastest, reaching about thirty-one percent error at a ratio of twenty. Population weighting reaches about fourteen percent and dasymetric about nine percent at the same ratio. The message is that no method rescues a transfer from a very coarse source. No method rescues a transfer from a very coarse source 0% 20% 40% median error area weighting population weighting dasymetric 1× 10× 20× source zone area ÷ target zone area Above about 5×, publish on the source geography instead and say why

The record that makes a transfer auditable is short, and it is worth being explicit about its four parts:

What the Crosswalk Record Has to Contain The four items that make an interpolation reproducible: the crosswalk fractions with stable source and target identifiers, the weighting layer identity and vintage, the method name with its parameters, and the conservation residual. Each is small, all four together are the only way to reproduce, revise or reverse the transfer, and none can be reconstructed from the interpolated values alone. Four items, none reconstructable from the output Crosswalk src_id, tgt_id, frac one row per pair the analysis itself Weighting layer identity + vintage blocks, 2020 decides the answer Method name + parameters population-weighted what to compare against Residual conservation check 4.2e−12 drift detector Persist all four with the output, keyed by run signature regenerating the crosswalk later from re-downloaded boundaries will not reproduce it

Compliance & Audit Controls

  • Persist the crosswalk, not just the result. The fractions are the analysis; regenerating them later from a changed boundary file will not reproduce the published numbers.
  • Record the weighting layer and its vintage in the run metadata beside the CRS, since it determines the answer as much as the method does.
  • Label interpolated columns as interpolated. A tract count that was measured and one that was distributed must not share a column name, because a downstream analyst cannot tell them apart and will treat both as observed.
  • Re-run disclosure review after interpolation. Splitting a count across targets creates new small cells that did not exist in the source, and the suppression decision has to be made on the published geography.
  • State the size ratio. If the source zones average more than about five times the target zones, say so next to the map; it is the single number that most constrains how much the result should be trusted.

When Not to Interpolate

The strongest recommendation in this whole topic is negative, and it is worth stating on its own because it is the option most often skipped. In several common situations the correct action is to publish on the geography the data arrived on and explain why.

When the size ratio is extreme. Above roughly five source zones per target zone, the transfer distributes an assumption rather than moving information, and the resulting fine-grained map has the visual authority of measurement with none of its content. Publishing the coarse map and saying “this is the resolution the data supports” is both more honest and, in practice, more persuasive to a reviewer than a smooth surface nobody can defend.

When the target analysis is a cluster test. Interpolated values within one source zone are perfectly correlated by construction, so any neighbour-based statistic computed on them will find structure that the interpolation put there. If the analysis plan ends in a hotspot map, either obtain the data at the analysis geography or change the analysis.

When the counts are small. A count of three distributed across five targets produces five fractions, none of which is a case, and the arithmetic that follows — rates, thresholds, suppression decisions — treats those fractions as though they were observations. Aggregate up instead of splitting down.

When the point data exists somewhere. This is the most frequent case and the easiest to miss, because the aggregated file is the one that arrived and the point file is one department away. Asking costs an email; interpolating costs an assumption that then travels through every downstream product.

A useful discipline is to record, in the analysis plan, the condition under which interpolation would be abandoned — the size ratio, the count floor, the validation misplacement share — before any transfer is run. Deciding afterwards that the result is good enough is not a decision; it is a description of what happened.

Implementation Checklist

FAQ

Can I interpolate a rate directly if I have no denominator? Only as a last resort, and it must be labelled. Interpolating a rate treats every part of the source zone as equally populated, which is the assumption population weighting exists to avoid. If the denominator is genuinely unavailable, an area-weighted rate transfer is defensible for display and not for analysis.

Is dasymetric interpolation always better? It is better when the ancillary layer is current and correctly classified, and worse when it is not, because a misclassified land-use polygon moves population confidently in the wrong direction. Its advantage over population weighting is also small when the source zones are already small.

How do I interpolate between two geographies that both changed? Do it in two steps through a common fine layer — usually blocks — rather than directly. A direct crosswalk between two moved geographies compounds two sets of assumptions in a way that cannot be diagnosed afterwards.

What if the source zones overlap? Then they are not a partition and no crosswalk is valid. Resolve the overlap first, either by taking a hierarchy of precedence or by treating the overlap as its own source zone with an apportioned count.

Recording a Transfer So Somebody Can Undo It

Interpolation is one of the few preparation steps that cannot be reversed from its output alone, so the record has to carry enough to reconstruct it. Four items suffice, and all four are cheap to persist at the moment the transfer runs: the crosswalk fractions with source and target identifiers, the weighting layer’s identity and vintage, the method name and any parameters it took, and the conservation check’s residual. With those, a later analyst can reproduce the numbers exactly, apply a different method to the same sources, or roll the transfer back to publish on the source geography. Without them, the interpolated values are a dead end, and the most common consequence is that the whole transfer is redone from scratch against a boundary file that has since changed.