Spatial Epi Public Health GIS

Dasymetric Refinement with Land-Use Ancillary Data

Areal-weighted interpolation assumes population is spread uniformly across a polygon, so a rate denominator built that way is wrong wherever a unit contains lakes, forest, or industrial land — and any suppression decision resting on that denominator is wrong with it. This guide, part of dasymetric suppression and small-count disclosure control, redistributes population onto inhabited land using a land-use mask and re-aggregates it to target zones with a mass-preserving equal-area overlay.

Problem Context & Constraints

Areal weighting reallocates a source count to a target zone in proportion to the overlap area: if 30% of a source tract’s area falls inside a target hexagon, the tract contributes 30% of its population. The method is simple and mass-preserving, but it encodes one assumption that fails constantly in health geography — that people are distributed uniformly across every square metre of the source polygon. A county that is 70% national forest does not have 70% of its residents living among the trees. The uniform assumption inflates the denominator of any target zone that captures the empty land and deflates the zone that captures the town, so the mapped rate is understated where people actually live.

For a source zone ss and target zone tt, plain areal weighting estimates the target population as

P^t=∑sPs As∩tAs,\hat{P}_t = \sum_s P_s \, \frac{A_{s \cap t}}{A_s},

where As∩tA_{s \cap t} is the intersection area and AsA_s the source area. Dasymetric refinement replaces raw area with inhabited area drawn from an ancillary layer:

P^t=∑sPs ∑cwc As∩tc∑cwc Asc,\hat{P}_t = \sum_s P_s \, \frac{\sum_{c} w_c\, A^{c}_{s \cap t}}{\sum_{c} w_c\, A^{c}_{s}},

where AcA^{c} is the area of land-use class cc and wcw_c is its residential weight — 00 for water and 11 for residential in the binary case, or graded weights (residential >> commercial >> agriculture >> 0 for water) in the weighted case. The denominator normalizes per source zone, which is what preserves total mass: every source resident is redistributed, never created or lost.

Two constraints shape the implementation:

  • Area must be measured on an equal-area CRS. Weights are ratios of areas, and areas computed in degrees distort with latitude, so the intersection must be projected before any .area call. This is the same projection discipline enforced across coordinate reference systems for public health.
  • Mass must be conserved exactly. The redistributed grand total must equal the source grand total to floating-point tolerance, or the denominator you hand to a rate map is silently wrong. This is a hard validation gate, not a nicety.
Dasymetric Reallocation Across a Land-Use Mask Three panels left to right. The left panel is a source zone with a uniform population count. The middle panel overlays a land-use mask that marks part of the zone as inhabited residential land and the rest as uninhabited water and forest. The right panel shows two target zones receiving population only in proportion to the inhabited land they capture, so the target that overlaps the uninhabited area receives little population while the target over the town receives most of it. Total population is unchanged across the panels. Population follows inhabited land, not polygon area — total mass is preserved 1 · Source zone P = 1000 assumed uniform 2 · Land-use mask inhabited w = 1 water forest w = 0 3 · Target zones town: 880 most inhabited land forest: 120 little inhabited land 880 + 120 = 1000 — the ancillary mask redistributes mass, it never changes the total Areal weighting would split 1000 by area alone, over-filling the forest target

Prerequisites

Pin the stack so the redistributed denominator is reproducible:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x and pyproj 3.6.x)
  • numpy 1.26.4, pandas 2.2.2

Input state assumed below:

  • Source zones: a GeoDataFrame of population polygons (census tracts) with a stable geoid and a population column.
  • Target zones: the geography you want denominators for (hex grid, service areas), with a stable target_id.
  • Ancillary land-use layer: inhabited-land polygons. In practice these come from vectorizing NLCD-style raster classes (residential/developed classes kept, water/wetland/forest dropped) or from a building-footprint layer. Each polygon carries a weight — 1 for a binary inhabited mask, or graded residential weights. How such a raster becomes an analysis-ready vector layer is covered in spatial data types and formats.
  • CRS: everything is projected to a single equal-area CRS (EPSG:5070 CONUS Albers) via pyproj before any area computation.

Step-by-Step Solution

The routine intersects source zones with the weighted land-use mask to get each source’s inhabited “mass units,” intersects again with the target zones, and apportions each source population by its share of inhabited mass falling in each target. The two-way overlay is what lets population cross from source geography to target geography through the mask.

# Mass-preserving dasymetric reallocation onto inhabited land.
# Pinned: python==3.11, geopandas==0.14.4, numpy==1.26.4, pandas==2.2.2
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("dasymetric")

EQUAL_AREA_CRS = "EPSG:5070"  # CONUS Albers — area ratios must be equal-area

def dasymetric_reallocate(source: gpd.GeoDataFrame,
                          target: gpd.GeoDataFrame,
                          landuse: gpd.GeoDataFrame,
                          *, pop_col="population") -> pd.DataFrame:
    """Redistribute source population to target zones via a weighted land-use mask.

    Returns a DataFrame of target_id -> dasymetric_pop, preserving total mass.
    """
    src = source.to_crs(EQUAL_AREA_CRS)[["geoid", pop_col, "geometry"]]
    tgt = target.to_crs(EQUAL_AREA_CRS)[["target_id", "geometry"]]
    lu = landuse.to_crs(EQUAL_AREA_CRS)[["weight", "geometry"]]

    # 1. inhabited mass = weighted area of (source ∩ land-use), per source zone
    src_lu = gpd.overlay(src, lu, how="intersection", keep_geom_type=True)
    src_lu["mass"] = src_lu.geometry.area * src_lu["weight"]
    src_mass = src_lu.groupby("geoid")["mass"].sum().rename("src_mass")

    # 2. inhabited mass in each (source ∩ target ∩ land-use) piece
    #    overlay the source-landuse pieces onto the target zones
    pieces = gpd.overlay(src_lu, tgt, how="intersection", keep_geom_type=True)
    pieces["piece_mass"] = pieces.geometry.area * pieces["weight"]

    # 3. each source's population follows its share of inhabited mass
    pieces = pieces.merge(src_mass, on="geoid", how="left")
    pieces = pieces.merge(src[["geoid", pop_col]], on="geoid", how="left")
    share = np.divide(pieces["piece_mass"], pieces["src_mass"],
                      out=np.zeros(len(pieces)),
                      where=pieces["src_mass"] > 0)
    pieces["alloc_pop"] = pieces[pop_col] * share

    out = (pieces.groupby("target_id")["alloc_pop"].sum()
           .rename("dasymetric_pop").reset_index())
    log.info("reallocated %d source zones onto %d target zones",
             src["geoid"].nunique(), out["target_id"].nunique())
    return out

For a binary mask, set every weight to 1 and drop non-residential polygons before calling. For a weighted scheme, assign graded weights per class — the normalization by src_mass guarantees mass preservation regardless of the weight scale, so the weights only have to be relatively correct, not absolute densities.

The choice between a binary and a weighted mask is a defensibility trade-off, not a purely technical one. A binary inhabited/uninhabited mask makes exactly one claim — that population is uniform within inhabited land and zero outside it — and that claim is easy to justify from a single land-cover reclassification. A weighted scheme (residential >> mixed-use >> commercial >> agriculture >> 0) is more realistic where developed density varies, but every weight is now an assumption a reviewer can challenge, and the weights propagate into every rate the denominator supports. The conservative default is binary; escalate to weighted only when you can point to an empirical density gradient for the classes involved and are prepared to log the exact mapping. When NLCD-style classes are the source, a common binary rule keeps the developed classes (open-space through high-intensity developed) and drops water, wetland, forest, and cultivated land — but “developed open space” is parks and golf courses as often as it is low-density housing, so treat that class with caution or exclude it.

Vectorizing the raster before overlay, rather than sampling the raster per polygon, is a deliberate choice here: an exact polygon-on-polygon intersection preserves mass to floating-point tolerance, whereas pixel counting introduces edge-quantization error that the mass-preservation assertion below would flag as a leak. The trade is cost — the double overlay is the expensive step — which is why target geographies are best precomputed once and cached rather than rebuilt per run.

Validation & Edge Cases

The single most important check is mass preservation: the redistributed grand total must equal the source grand total. Assert it, so a leak halts the pipeline.

def assert_mass_preserved(source, allocated, *, pop_col="population",
                          rtol=1e-6):
    src_total = float(source[pop_col].sum())
    alloc_total = float(allocated["dasymetric_pop"].sum())
    log.info("mass check: source=%.3f allocated=%.3f", src_total, alloc_total)
    if not np.isclose(src_total, alloc_total, rtol=rtol):
        raise AssertionError(
            f"mass not preserved: {src_total:.3f} != {alloc_total:.3f}; "
            "a source zone with no inhabited land likely lost its population")

Two edge cases break mass preservation if unhandled:

1. Source zones with no inhabited land. If a source zone’s geometry never intersects the land-use mask — a wholly protected watershed that nonetheless has a recorded population, or a mask with a coverage gap — its src_mass is zero, the share divide yields zero everywhere, and that zone’s population silently vanishes. The mass-preservation assertion catches it. The fix is a documented fallback: for zero-inhabited-land source zones, fall back to plain areal weighting so their residents are still placed somewhere, and log the fallback explicitly.

zero_land = set(source["geoid"]) - set(src_mass[src_mass > 0].index)
if zero_land:
    log.warning("%d source zones have no inhabited land; areal fallback applied: %s",
                len(zero_land), sorted(zero_land)[:10])

2. Mask coverage gaps at the study edge. A target zone straddling the study-area boundary captures inhabited land that has no source population on the far side, so its denominator is understated. Pad the source and mask layers one zone beyond the reporting boundary and clip the target output back at the end, so edge zones are apportioned against complete inhabited land.

A third, quieter failure is a CRS mismatch that survives silently: if the land-use layer arrives in a geographic CRS and the reprojection is skipped, .area returns square degrees and the weights are still ratios, so mass is preserved but the distribution is wrong. Guard it by asserting the working CRS is projected and equal-area before the first overlay:

assert source.to_crs(EQUAL_AREA_CRS).crs.is_projected, "reproject to equal-area first"

This is the dangerous class of bug precisely because the mass-preservation check passes: ratios of square degrees still sum to one per source zone, so the grand total is conserved while individual target denominators are quietly skewed by the latitude-dependent distortion. Only an equal-area assertion catches it, which is why the guard is an assertion and not a comment. A related sanity check worth logging is the implied density per target zone (dasymetric_pop / inhabited_area): a target whose implied density is orders of magnitude off its neighbours usually points to a mask alignment error or a stray high-weight polygon, not a real settlement pattern.

One more edge case deserves explicit handling: a target zone with inhabited land but no source population, the mirror of the zero-inhabited-land case. This happens when the source population layer has a coverage gap the mask does not — a newly incorporated area, or a boundary vintage mismatch. Such a target simply receives zero allocated population, which is correct given the inputs but may not be true; log target zones that end with dasymetric_pop == 0 despite non-trivial inhabited area so the gap is visible rather than silently mapped as an empty denominator that later divides a rate to NaN.

Compliance Notes

A dasymetric denominator is a modeling choice that flows into every rate and every suppression decision downstream, so it must be reconstructable. Record three things in the run’s provenance. Log the ancillary source identity and vintage — which land-use product, which year, which classes were treated as inhabited — because a denominator built on 2011 land cover applied to 2020 population is a defect a reviewer must be able to see. Log the weight scheme: the exact class-to-weight mapping (binary, or the graded weights) belongs in the metadata, since it determines the redistribution. And log the fallback events: every source zone that reverted to areal weighting for lack of inhabited land is a place where the refinement did not apply, and honest reporting of those zones is what separates a defensible denominator from a black box. Bind these to a configuration hash so the reallocation, and therefore the rates and suppression that consume it, can be reproduced exactly. The downstream small-count cell suppression in choropleth rate maps step, and the rate-stability flags that ride on this denominator, inherit whatever error the ancillary layer carries, so its documentation is not optional.