⚙️ How-To Census & demographics Interpolation & surfaces Raster & rasterio Census & Demographics · step 14 of 20 Updated

How to Redistribute Population with Dasymetric Mapping

Problem statement

Population is published for census units, but questions are asked about other zones: a flood extent, a school catchment, a grid, a service area. The usual answer is area weighting, which splits a tract's people across the zones it overlaps in proportion to area. That assumes people are spread evenly across the tract, and they are not.

This guide tests the assumption where the answer is known. Delaware's 262 census tracts were redistributed to its 706 block groups, whose 2020 counts come from the census blocks. Area weighting put 202,169 people (20.4%) in the wrong block group, and the median block group was off by 34.6%.

Dasymetric mapping uses ancillary data to say where people can live within each source zone. A land-cover mask cut the error to 15.9%. Building footprints cut it to 11.0%. Removing water alone, the most common "improvement", took it from 20.4% to 20.1%.

Quick answer

With a land-cover raster, tobler does the whole thing in one call:

from tobler.dasymetric import masked_area_interpolate

estimate = masked_area_interpolate(
    tracts.to_crs(5070),                 # source zones with the counts
    block_groups.to_crs(5070),           # target zones
    raster="worldcover_30m.tif",         # ESA WorldCover, resampled to 30 m
    pixel_values=[50],                   # class 50 = built-up
    extensive_variables=["POP20"],
    nodata=0,
)

Measured on Delaware: 10.8 s, total preserved at 989,948, and 165,096 people misallocated (16.7%) against 202,169 for plain area weighting.

With building footprints instead of a raster, the dasymetric() function in Example 1 took 0.9 s and misallocated 108,530 people (11.0%).

Bar chart of population misallocated: area weighting 20.4%, land only 20.1%, tobler mask 16.7%, WorldCover built-up 15.9%, building footprint 11.5% and building count 11.0%.
The ancillary data that knows where homes are halves the error of plain area weighting.

Step-by-step solution

1. Redistribute counts, not rates

Dasymetric mapping splits an extensive quantity — people, households, housing units — across space. Medians, percentages and densities do not split. Redistribute the counts behind them and recompute on the target zones.

2. Build a test you can score

Pick a pair of geographies where the target totals are known. Census blocks sum exactly to block groups and tracts, so you can dissolve blocks into both and pretend you only had the tracts:

import geopandas as gpd

blocks = gpd.read_file("tl_2020_10_tabblock20.zip", columns=["GEOID20", "POP20"])
tracts = blocks.dissolve(by=blocks.GEOID20.str[:11], aggfunc={"POP20": "sum"})
groups = blocks.dissolve(by=blocks.GEOID20.str[:12], aggfunc={"POP20": "sum"})
truth = groups["POP20"]

The TIGER/Line 2020 block file already carries POP20, so no API call is needed. Dissolving Delaware's 20,198 blocks took 2.7 s.

3. Measure the baseline first 

from tobler.area_weighted import area_interpolate

aw = area_interpolate(tracts.to_crs(5070), groups.to_crs(5070),
                      extensive_variables=["POP20"])

0.45 s, and 20.4% of people in the wrong block group. That is the number any ancillary dataset has to beat. Without it you cannot tell whether a mask helped.

4. Choose the ancillary data

Two open options cover most of the world:

  • Land cover. ESA WorldCover 2021 is a 10 m global classification, published as cloud-optimised GeoTIFF tiles on AWS. Class 50 is built-up. Two tiles cover Delaware.
  • Building footprints. Microsoft's US Building Footprints has 357,534 buildings in Delaware. The state file is a 13.5 MB zip and read in 2.3 s.

Both are more recent than the 2020 census, so a few buildings will postdate the counts. That error is small compared with the one area weighting makes.

5. Mask with land cover in tobler

masked_area_interpolate turns the matching raster pixels into polygons, clips each source zone to them, then area-weights the clipped zones. Vectorising 10 m pixels is expensive: the Delaware window held 6.4 million built-up pixels. Resample to 30 m with a mode filter first:

import rasterio
from rasterio.enums import Resampling

with rasterio.open("worldcover_2021_delaware.tif") as src:
    h, w = src.height // 3, src.width // 3
    data = src.read(1, out_shape=(h, w), resampling=Resampling.mode)
    transform = src.transform * src.transform.scale(src.width / w, src.height / h)
    profile = src.profile | {"height": h, "width": w, "transform": transform, "nodata": 0}
with rasterio.open("worldcover_30m.tif", "w", **profile) as dst:
    dst.write(data, 1)

The 30 m raster had 693,620 built-up pixels, and the tobler call took 10.8 s with four jobs. Pass nodata=0 for WorldCover, because tobler's default is 255.

6. Weight by buildings

Buildings are closer to homes than a built-up pixel is: a built-up pixel can be a car park or a motorway junction. Use each building's representative point, count them in each source × target intersection, and share each source's people by those counts (Example 1):

ancillary weight                  people misallocated   median |error|
none (area weighting)                202,169  (20.4%)            34.6%
land only (water removed)            199,102  (20.1%)            33.1%
WorldCover built-up, 10 m            157,841  (15.9%)            25.9%
tobler mask, WorldCover 30 m         165,096  (16.7%)            27.4%
building count                       108,530  (11.0%)            16.7%
building footprint area              113,377  (11.5%)            18.6%
buildings of 40–400 m² only          121,355  (12.3%)            18.2%

7. Test a target that does not nest

Block groups sit inside tracts, which flatters every method. A 2 km grid does not nest. Scored against block populations assigned to cells by their representative points:

2 km grid, 3,773 cells     misallocated     people put in cells that have nobody
area weighting          226,435 (22.9%)                               18,463
building count          117,920 (11.9%)                                3,041

The last column is what a map shows: area weighting paints population into fields, marsh and forest.

8. Keep the totals and a fallback

Every method above preserved the state total of 989,948. A source zone with no ancillary features at all, such as a tract that is all water or park, must fall back to area weighting. Otherwise its people vanish. Two Delaware block groups had no buildings; both were empty in the census too.

Table comparing area weighting and building-count dasymetric mapping onto a 2 km grid: people misallocated and people placed in cells with no census population.
On a grid, area weighting paints people into fields and forest.

Code examples

Example 1 — point-weighted dasymetric redistribution for any target zones 

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


def dasymetric(source, target, value, ancillary_points, weight=None, crs=5070):
    """Split source totals across target zones in proportion to ancillary points.

    Works for targets that do not nest in the sources: the allocation happens on
    the source x target intersections. Where a source has no ancillary points at
    all, it falls back to area weighting so nobody is dropped.
    """
    src = source[[value, "geometry"]].to_crs(crs).reset_index(names="src_id")
    tgt = target[["geometry"]].to_crs(crs).reset_index(names="tgt_id")
    pieces = gpd.overlay(src, tgt, how="intersection", keep_geom_type=True)
    pieces["area"] = pieces.area

    pts = ancillary_points.to_crs(crs)
    pts = pts.assign(_w=pts[weight] if weight else 1.0)[["_w", "geometry"]]
    hits = gpd.sjoin(pts, pieces[["geometry"]], predicate="within", how="inner")
    pieces["anc"] = hits.groupby("index_right")["_w"].sum().reindex(pieces.index).fillna(0)

    by_source = pieces.groupby("src_id")[["anc", "area"]].transform("sum")
    share = np.where(by_source["anc"] > 0, pieces["anc"] / by_source["anc"].where(by_source["anc"] > 0),
                     pieces["area"] / by_source["area"])
    pieces["estimate"] = pieces[value] * share
    out = pieces.groupby("tgt_id")["estimate"].sum()
    return out.reindex(tgt["tgt_id"]).fillna(0).set_axis(target.index)
buildings = gpd.read_file("/vsizip/Delaware.geojson.zip").to_crs(5070)
points = gpd.GeoDataFrame({"footprint": buildings.area},
                          geometry=buildings.representative_point(), crs=5070)
by_count = dasymetric(tracts, groups, "POP20", points)
by_area = dasymetric(tracts, groups, "POP20", points, weight="footprint")

Example 2 — counting land-cover classes per zone without vectorising pixels 

import rasterio
from rasterio.features import rasterize


def raster_class_counts(zones, raster_path, classes):
    """Count pixels of the given classes inside each zone (zones in any CRS)."""
    with rasterio.open(raster_path) as src:
        data = src.read(1)
        shapes = zones.to_crs(src.crs).geometry
        ids = rasterize(((geom, i + 1) for i, geom in enumerate(shapes)),
                        out_shape=data.shape, transform=src.transform,
                        fill=0, dtype="uint32")
    mask = np.isin(data, classes)
    counts = np.bincount(ids[mask], minlength=len(zones) + 1)[1:]
    return pd.Series(counts, index=zones.index)

Burning the 706 block groups into the 16,920 × 9,120 pixel window and counting built-up pixels took 1.1 s, against 10.8 s for tobler's vectorising route. For zones that nest in the sources, the counts become weights directly. For non-nesting zones, pass source × target intersections as zones. Each rasterised zone keeps the pixels whose centres fall inside it, so very small pieces can receive none.

Example 3 — scoring a redistribution 

def score(estimate, truth):
    """Error of a redistribution against known target totals."""
    error = estimate.reindex(truth.index).fillna(0) - truth
    pct = error.abs() / truth.where(truth > 0) * 100
    return {
        "total": round(float(estimate.sum())),
        "misallocated": round(float(error.abs().sum() / 2)),
        "misallocated_pct": round(float(error.abs().sum() / 2 / truth.sum() * 100), 1),
        "median_abs_pct": round(float(pct.median()), 1),
    }
area-weighted        {'total': 989948, 'misallocated': 202169, 'misallocated_pct': 20.4, 'median_abs_pct': 34.6}
buildings (count)    {'total': 989948, 'misallocated': 108530, 'misallocated_pct': 11.0, 'median_abs_pct': 16.7}

"Misallocated" is half the sum of absolute errors. Each misplaced person appears once as a surplus and once as a deficit.

Explanation

Why area weighting fails inside a single tract

A tract is drawn to hold about 4,000 people, not to be uniform. On its edge, a suburban tract may include a subdivision, a golf course and farmland. Area weighting gives the farmland its share of the subdivision's residents.

The worst Delaware block group, 100030145021, has 2,312 residents. Area weighting gave it 5,977 because it holds a large share of its tract's land. The building-weighted estimate was 4,059: still wrong, but less than half as wrong.

Why removing water barely helped

Water masking is the textbook first step, and here it moved the error from 20.4% to 20.1%. Delaware's water is concentrated in a few coastal and bay tracts. Inside a typical tract, the people are unevenly spread over dry land. A mask only helps if it removes the area where people are not, and on land that means knowing where the buildings are.

Why buildings beat land cover

The built-up class marks impervious surface, and impervious surface includes car parks, industrial estates and roads. It improved on area weighting (15.9%) but still spread residents over commercial land.

Building counts come closer to dwellings. Counting buildings (11.0%) beat weighting by footprint area (11.5%). A warehouse's large footprint would otherwise attract the people who live in the houses nearby. Keeping only buildings of 40–400 m² made things worse (12.3%). The filter discards terraces and apartment blocks mapped as single large polygons, which are exactly where people are densest.

Why the grid result matters more than the block group result

Targets that nest inside the sources limit how wrong any method can be, because each source's people stay inside it. A grid cuts across tracts, and there area weighting put 18,463 people in cells where the census counted none. On a population grid, those phantom residents are the first thing a reader notices.

What the remaining 11% is

Buildings are a proxy. Garages, sheds, barns, shops and schools are counted as if they were homes. Vacant homes are counted as occupied. The footprints are from years after the 2020 census. Better proxies include residential address points, building type attributes and housing unit counts, and each of those would reduce the error further. Measure before assuming any of them does.

Bar chart for block group 100030145021: 2,312 people counted, 5,977 by area weighting and 4,059 by building count.
Buildings halved the error here, and a proxy that is not a home still leaves some behind.

Edge cases or notes

  • Project before measuring area. tobler warns "Geometry is in a geographic CRS" if you pass EPSG:4269 zones; the result here moved by at most 1.7 people, but do not rely on that elsewhere.
  • tobler's nodata defaults to 255. WorldCover uses 0; pass it explicitly.
  • Resample large rasters before masked_area_interpolate. It vectorises every matching pixel before overlaying.
  • WorldCover is in EPSG:4326. Pixel area varies with latitude by about 2% across a state the size of Delaware; for a continent, weight by pixel area.
  • Use representative points, not centroids, for footprints. A centroid of an L-shaped building can fall in the neighbouring zone.
  • 8,539 Delaware footprints fell outside every block group, mostly beyond the census coastline; they cannot receive people and should not.
  • Fallback is required. A source zone with no ancillary features must keep its people somewhere.
  • Intensive variables do not redistribute. tobler's intensive_variables averages them; it does not make a median meaningful on a new zone.

FAQ

What is dasymetric mapping?

Redistributing a count from source zones to other zones using ancillary data about where the count can occur, such as buildings or built-up land. It is area weighting restricted to the places people actually live.

How much better than area weighting is it?

On Delaware tracts redistributed to block groups, area weighting misallocated 20.4% of the population, a WorldCover built-up mask 15.9% and building counts 11.0%. On a 2 km grid, buildings halved the error, from 22.9% to 11.9%.

Which Python function does it?

tobler.dasymetric.masked_area_interpolate takes a raster and the pixel values to keep. For vector ancillary data such as building footprints, weight the source × target intersections yourself, as in Example 1.

Is removing water enough?

Not in this test: it improved the error from 20.4% to 20.1%. The unevenness that matters is on land, between built and unbuilt areas.

Should I weight by building count or footprint area?

Count did slightly better here (11.0% against 11.5%), because large non-residential footprints attract people under area weighting. Filtering to small buildings made it worse.

How do I know it worked on my data?

Score it where the truth is known. Dissolve census blocks into two levels, redistribute from the upper to the lower, and compare against the block totals before trusting the method on zones where you cannot check.