Area-Weighted Interpolation for Partial Block Group Overlap

This page solves one exact task: when a store catchment polygon covers only part of a census block group, apportioning that block group’s population to the catchment by the fraction of its area that falls inside, rather than counting all of it or none of it.

A drive-time catchment almost never lands on tidy block-group boundaries. Its edge slices through dozens of block groups, and a naive containment join forces a binary choice for each one — include the whole population or discard it. Both answers are wrong at the boundary, and the errors do not cancel; they compound into a trade-area total that can miss by 20% or more on a compact urban catchment. Areal interpolation replaces that binary with a fraction, and turns a coarse count into a defensible estimate.

Prerequisites #

Before running this task you need:

• Python packages: geopandas, shapely, and pandas. Install with pip install geopandas shapely pandas.
• A block-group layer with a count variable. Population, households, or any extensive (additive) quantity, joined to block-group geometry. Retrieving these is the upstream task in how to join ACS 5-year estimates to custom trade area polygons.
• A catchment polygon. A single trade-area geometry — typically a drive-time contour — in a known CRS.
• Valid input geometry. Overlay is unforgiving of self-intersections and unclosed rings; repair invalid geometry before intersecting, or the overlay silently drops features. Cleaning boundary artifacts is covered in fixing sliver polygons in spatial join operations.

The Areal-Weighting Assumption and Its Limits #

Areal interpolation rests on one assumption: that the target quantity is distributed uniformly across each source block group. Under that assumption, the share of a block group’s population belonging to the catchment equals the share of its area that the catchment covers. For block group $i$ with area $A_i$ and intersection with catchment $T$ of area $A_{i \cap T}$, the area weight is:

$$w_i = \frac{A_{i \cap T}}{A_i}, \qquad 0 \le w_i \le 1$$

The interpolated total over all intersecting block groups is the population-weighted sum:

$$\hat{P}_T = \sum_i w_i \, P_i$$

The uniform-density assumption is the method’s strength and its ceiling. Block groups are drawn to hold roughly equal population, so within a dense neighborhood the assumption is reasonable. It breaks where population is clustered inside a large, mostly-empty polygon — a block group spanning a reservoir, an industrial park, or a single apartment tower surrounded by parkland. There, half the area may hold none of the people, and pure area weighting misallocates them. When that error matters, dasymetric refinement (masking uninhabitable land, or weighting by a finer surface) improves on plain areal weighting, but adds inputs and complexity; area weighting is the correct, auditable default and the right starting point.

Configuration and Execution Parameters #

Why the CRS Is Not Optional Here #

Area weighting divides one area by another, and area is meaningless in a geographic (degree-based) CRS. In EPSG:4326 a “square degree” covers a wildly different amount of ground in Miami than in Minneapolis, so a weight computed there is distorted by latitude before you even start. Every area in this task — the block-group area $A_i$ and the intersection area $A_{i \cap T}$ — is computed in an equal-area projection so the ratio is honest. The code asserts the projection before any .area call.

Annotated Implementation #

The function reprojects both layers to an equal-area CRS, repairs geometry, records each block group’s full area, overlays the intersection to get the overlap area, computes the clamped weight, and returns both the apportioned per-block-group parts and the interpolated total.

import geopandas as gpd
import pandas as pd
from shapely.validation import make_valid
from pyproj import CRS

EQUAL_AREA_CRS = CRS.from_epsg(5070)   # CONUS Albers; areas are valid here
SLIVER_MIN_M2 = 1.0                    # overlap fragments below this are noise


def area_weighted_population(
    block_groups: gpd.GeoDataFrame,
    catchment: gpd.GeoDataFrame,
    value_col: str = "population",
) -> tuple[gpd.GeoDataFrame, float]:
    """
    Apportion an extensive block-group count to a catchment by fractional
    overlap area (areal interpolation under a uniform-density assumption).

    Returns the per-block-group apportioned parts and the interpolated total.
    """
    # 1. CRS discipline: area ratios are only meaningful in an equal-area CRS.
    assert block_groups.crs is not None, "block groups carry no CRS"
    assert catchment.crs is not None, "catchment carries no CRS"
    bg = block_groups.to_crs(EQUAL_AREA_CRS).copy()
    ta = catchment.to_crs(EQUAL_AREA_CRS).copy()

    # 2. Repair geometry so the overlay does not silently drop features.
    bg["geometry"] = bg.geometry.apply(make_valid)
    ta["geometry"] = ta.geometry.apply(make_valid)

    # 3. Record each block group's FULL area BEFORE clipping (the denominator).
    bg["bg_area"] = bg.geometry.area
    bg = bg[bg["bg_area"] > 0]

    # 4. Intersection overlay: keep only the polygonal overlap of the two layers.
    overlap = gpd.overlay(bg, ta, how="intersection", keep_geom_type=True)
    if overlap.empty:
        return overlap, 0.0

    # 5. Overlap area is the numerator; drop numerical slivers.
    overlap["overlap_area"] = overlap.geometry.area
    overlap = overlap[overlap["overlap_area"] >= SLIVER_MIN_M2]

    # 6. Area weight in [0, 1], then apportion the extensive count.
    overlap["weight"] = (overlap["overlap_area"] / overlap["bg_area"]).clip(0.0, 1.0)
    overlap["apportioned"] = overlap["weight"] * overlap[value_col]

    interpolated_total = float(overlap["apportioned"].sum())
    return overlap, interpolated_total

A single call yields the total plus an auditable per-block-group breakdown:

parts, total = area_weighted_population(bg_with_pop, drive_time_catchment, "population")
print(f"interpolated catchment population: {total:,.0f}")
print(parts[["GEOID", "weight", "population", "apportioned"]].sort_values("weight").head())

Three details carry the correctness of this function. First, bg_area is captured before the overlay, because after the intersection the geometry is already clipped and .area would return the overlap, forcing every weight to 1.0. Second, keep_geom_type=True prevents an edge-touching intersection from yielding a LineString whose zero area corrupts the ratio. Third, the weight is clamped to [0, 1]; floating-point overlay can produce an overlap area a hair larger than the parent, and an unclamped weight above 1.0 would over-count.

Failure Modes and Debugging #

Slivers deserve emphasis: an edge-aligned boundary between the catchment and a block group can produce hairline overlap polygons of a few square meters, which contribute negligible population but clutter the output and can skew a conservation check. The area floor discards them, and the broader treatment of these artifacts lives in fixing sliver polygons in spatial join operations.

Verification #

Confirm the interpolation before it feeds a score:

1. Weights in range: assert parts["weight"].between(0, 1).all() — any weight outside [0, 1] signals a CRS or denominator bug.
2. Conservation check: a block group fully inside the catchment must contribute its entire population. Filter to fully-contained block groups and assert weight ≈ 1.0; a block group fully outside must not appear at all.
3. Bounded total: the interpolated total must be less than or equal to the raw sum of every touched block group and greater than or equal to the sum of only the fully-contained ones. Outside that band, the weighting is wrong.
4. Overlap-area sanity: parts["overlap_area"].sum() should approximate the catchment’s own area (where the catchment lies over populated land), confirming the overlay covered the trade area without gaps or double coverage.

full = parts[parts["weight"] > 0.999]["population"].sum()          # fully inside
touched = parts["population"].sum()                                # any overlap
assert full <= total <= touched, "interpolated total outside conservation bounds"

catchment_area = drive_time_catchment.to_crs(EQUAL_AREA_CRS).geometry.area.sum()
assert parts["overlap_area"].sum() <= catchment_area * (1 + 1e-6), "overlap exceeds catchment"

Once populations are apportioned, they inherit the sampling uncertainty of the underlying ACS figures, so the same fractional weights should carry through to the margin-of-error propagation described in handling ACS margin of error in trade area estimates. The end-to-end mechanics of joining demographics to catchments — of which this apportionment is one refinement — are the subject of the demographic data integration and spatial joins reference.

Related #

• Fixing Sliver Polygons in Spatial Join Operations — clean the hairline overlap artifacts that pollute area weights.
• How to Join ACS 5-Year Estimates to Custom Trade Area Polygons — produce the population counts this method apportions.
• Demographic Data Integration & Spatial Joins — the full Census-to-catchment pipeline this task feeds.

← Back to Performing Point-in-Polygon Joins for Store Catchments