How to Interpolate to Polygons Instead of a Grid

Problem statement

Population by census tract, sales by postcode, cases by health district: the data arrives on one set of zones and the question is about another. Areal interpolation moves values between incompatible polygon sets.

The default method is area weighting: assume the value is spread evenly within each source zone and split it by overlap area. It is transparent, volume-preserving and usually wrong in a known direction — populations are not evenly spread.

The standard improvement is dasymetric weighting: use an ancillary variable such as building footprints to distribute the value more realistically. Measured on Leiden's 58,879 OpenStreetMap buildings, redistributing building counts from a 1 km grid to an offset 700 m grid and scoring both against the true counts (Example 4):

                  RMSE     MAE     bias    total
area weighting  351.47  258.33   -50.96   50,283
dasymetric      458.84  247.01  +102.03   58,851
truth                                     53,137

The dasymetric version was 31% worse on RMSE, while slightly better on MAE — 247.01 against 258.33. That is not a coding error; it is what happens when the ancillary variable is not proportional to the target: most cells come out a little better and the dense ones a lot worse.

Quick answer

import geopandas as gpd


def area_weighted(source, target, value_column, source_id="sid",
                  target_id="tid"):
    """Split each source value by the share of its area in each target zone."""
    pieces = gpd.overlay(
        source[[source_id, value_column, "geometry"]],
        target[[target_id, "geometry"]], how="intersection")

    source_area = source.set_index(source_id).geometry.area
    pieces["fraction"] = pieces.geometry.area / pieces[source_id].map(source_area)
    pieces["estimate"] = pieces[value_column] * pieces["fraction"]

    out = pieces.groupby(target_id)["estimate"].sum()
    print(f"  source total {source[value_column].sum():,.0f}, "
          f"redistributed {out.sum():,.0f}")
    return out.reindex(target[target_id], fill_value=0.0)

The two totals should match. A difference means source zones extend beyond the target zones — which is information, not a bug.

58,879 OSM buildings in Leiden, EPSG:4326; columns ['building', 'name', 'height', 'levels']

Extensive or intensive, area weights, then dasymetric

1. Decide whether the value is extensive or intensive

Extensive values — counts, totals, populations — split between zones and must be volume-preserving. A tract with 500 people contributing 40% of its area contributes 200 people.

Intensive values — densities, rates, means, percentages — do not split. They are area-weighted averages, and summing them is meaningless.

if intensive:
    estimate = (value * intersection_area).sum() / target_area
else:
    estimate = (value * area_fraction).sum()

Using the extensive formula on a rate produces a number that grows with zone size, which is a common and consequential error.

2. Area weighting first, and know its assumption

Area weighting assumes uniform density within each source zone. That assumption is always false and often harmless: measured here it produced a bias of −50.96 buildings per zone against a median zone count of 722.

Its virtue is transparency. There is one assumption, it is stated, and it is easy to reason about.

3. Dasymetric weighting only helps if the ancillary variable is proportional

The measurement is the point of this page. Weighting by building footprint area to redistribute building count made the RMSE worse — 458.84 against 351.47 — even though the mean absolute error improved slightly.

The reason is that footprint area is not proportional to count. One warehouse of 5,000 m² is one building; twenty terraced houses of 80 m² are twenty. The ancillary variable predicts floor space, and the target was a count.

Dasymetric mapping is powerful when the ancillary variable really is proportional — residential footprint area for population, cropland area for yield. It is worse than nothing when it is not.

4. Check the totals

area weighting  redistributed 50,283 against a source total of 58,879
dasymetric      redistributed 58,851

Area weighting lost 8,596 because the source grid extends beyond the target grid, so part of each edge cell's value falls outside every target zone. Dasymetric weighting normalised within each source zone, so it kept all but 28 by pushing edge value into the target zones that did exist; the 28 sit where a source cell's pieces inside the target grid hold no buildings and the function falls back to area weighting. The truth, 53,137, lies between the two.

Neither is wrong; they answer different questions. Decide which behaviour you want at the edges and state it.

5. Report the assumption with the result

An areal-interpolated value is a model output. The zone geometries, the method and the ancillary variable are all part of what the number means.

Leiden buildings reallocated between zone sets

The data used on this page

The buildings are every building in Leiden in OpenStreetMap — about 59,000 footprints with their building type — fetched with OSMnx through an Overpass mirror. The zones are built from them in Example 4: a 1 km source grid over the whole city, and a 700 m target grid offset from it and inset from the edges, both in the local UTM zone. Because every building is known, the true count in each target cell is known too, which is what makes the methods scoreable. The OpenStreetMap answers are read as the Overpass mirror gave them on 9 October 2026.

import io
import os
import requests
import zipfile

import geopandas as gpd
import osmnx as ox

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-areal-interpolation-python-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # overpass-api.de mirror
ox.settings.overpass_rate_limit = False
if not os.path.exists("buildings.parquet"):
    found = ox.features_from_place("Leiden, Netherlands", {"building": True}).reset_index()
    found = found[found.geom_type.isin(["Polygon", "MultiPolygon"])]
    keep = [c for c in ("osmid", "building", "name", "height", "building:levels") if c in found.columns]
    found[keep + ["geometry"]].rename(columns={"building:levels": "levels"}).to_parquet("buildings.parquet")
buildings = gpd.read_parquet("buildings.parquet")
print(f"{len(buildings):,} OSM buildings in Leiden, EPSG:{buildings.crs.to_epsg()}; "
      f"columns {[c for c in buildings.columns if c != 'geometry']}")

58,879 OSM buildings in Leiden, EPSG:4326; columns ['building', 'name', 'height', 'levels']

Example 1 — extensive and intensive together

import geopandas as gpd
import numpy as np


def areal_interpolate(source, target, extensive=(), intensive=(),
                      source_id="sid", target_id="tid"):
    """Redistribute both kinds of value, with the totals checked."""
    pieces = gpd.overlay(
        source[[source_id, *extensive, *intensive, "geometry"]],
        target[[target_id, "geometry"]], how="intersection")
    pieces["piece_area"] = pieces.geometry.area

    source_area = source.set_index(source_id).geometry.area
    target_area = target.set_index(target_id).geometry.area
    pieces["source_fraction"] = pieces["piece_area"] / \
        pieces[source_id].map(source_area)

    out = {}
    for column in extensive:
        pieces["_e"] = pieces[column] * pieces["source_fraction"]
        out[column] = pieces.groupby(target_id)["_e"].sum()
        total_in = source[column].sum()
        total_out = out[column].sum()
        print(f"  {column:16} {total_in:12,.0f} -> {total_out:12,.0f} "
              f"({total_out / total_in - 1:+.2%})")

    for column in intensive:
        pieces["_w"] = pieces[column] * pieces["piece_area"]
        numerator = pieces.groupby(target_id)["_w"].sum()
        denominator = pieces.groupby(target_id)["piece_area"].sum()
        out[column] = numerator / denominator
        print(f"  {column:16} area-weighted mean, "
              f"{out[column].min():.3f}..{out[column].max():.3f}")

    frame = gpd.GeoDataFrame(out).reindex(target[target_id])
    frame.index.name = target_id
    return frame

Printing the total in and out for every extensive column is the check that catches an overlay problem, a CRS mismatch or partial coverage. A loss of a few percent is usually edge effects; a loss of 40% is a bug.

Example 2 — dasymetric weighting, with a proportionality test

import geopandas as gpd
import numpy as np


def dasymetric(source, target, value_column, ancillary,
               source_id="sid", target_id="tid", test=True):
    """Redistribute by an ancillary variable, after checking it is proportional."""
    if test:
        joined = gpd.sjoin(ancillary[["geometry"]].assign(
            weight=ancillary.geometry.area),
            source[[source_id, value_column, "geometry"]],
            predicate="within")
        per_zone = joined.groupby(source_id)["weight"].sum()
        values = source.set_index(source_id)[value_column]
        common = per_zone.index.intersection(values.index)
        correlation = float(np.corrcoef(per_zone[common], values[common])[0, 1])
        print(f"  correlation between the ancillary total and {value_column}: "
              f"{correlation:+.3f}")
        if correlation < 0.7:
            print("  ! weak proportionality — dasymetric weighting may be "
                  "worse than plain area weighting")

    pieces = gpd.overlay(
        source[[source_id, value_column, "geometry"]],
        target[[target_id, "geometry"]], how="intersection").reset_index()

    within = gpd.overlay(ancillary[["geometry"]],
                         pieces[["index", "geometry"]], how="intersection")
    within["weight"] = within.geometry.area
    piece_weight = within.groupby("index")["weight"].sum()

    pieces["weight"] = pieces["index"].map(piece_weight).fillna(0.0)
    total = pieces.groupby(source_id)["weight"].transform("sum")
    pieces["fraction"] = np.where(
        total > 0, pieces["weight"] / total,
        pieces.geometry.area / pieces[source_id].map(
            source.set_index(source_id).geometry.area))
    pieces["estimate"] = pieces[value_column] * pieces["fraction"]

    out = pieces.groupby(target_id)["estimate"].sum()
    print(f"  source total {source[value_column].sum():,.0f}, "
          f"redistributed {out.sum():,.0f}")
    return out.reindex(target[target_id], fill_value=0.0)

Run on the measured case, the test passed — a correlation of +0.874 between footprint area and building count per source cell — and the method still lost on RMSE. A correlation across zones shows that the ancillary variable tracks the value between zones; dasymetric weighting relies on proportionality within them, which a zone-level correlation cannot see. Treat the test as necessary, not sufficient; validate as in Example 3.

The fallback to area weighting where a source zone contains no ancillary features is essential. Without it, those zones lose their value entirely.

Example 3 — validating against known truth

import numpy as np


def validate_areal(estimates, truth, labels=None):
    """Compare methods against a target-zone count you actually have."""
    for name, estimate in estimates.items():
        error = np.asarray(estimate) - np.asarray(truth)
        print(f"  {name:18} RMSE {np.sqrt((error ** 2).mean()):8.2f}  "
              f"MAE {np.abs(error).mean():8.2f}  "
              f"bias {error.mean():+8.2f}  "
              f"total {np.sum(estimate):10,.0f} vs {np.sum(truth):,.0f}")

    best = min(estimates, key=lambda k: np.sqrt(
        ((np.asarray(estimates[k]) - np.asarray(truth)) ** 2).mean()))
    print(f"  best on RMSE: {best}")
    return best

Where the target-zone values are known for one period — a census year, a survey — validating there and applying the chosen method elsewhere is far better than assuming the sophisticated method wins.

Example 4 — the measurement

import numpy as np
from shapely.geometry import box

buildings = gpd.read_parquet("buildings.parquet")
buildings = buildings.to_crs(buildings.estimate_utm_crs())


def grid(bounds, size, offset=0.0, inset=0.0, id_name="sid"):
    x0, y0, x1, y1 = bounds
    x0, y0, x1, y1 = x0 + inset + offset, y0 + inset + offset, x1 - inset, y1 - inset
    cells = [box(x, y, x + size, y + size)
             for x in np.arange(x0, x1, size) for y in np.arange(y0, y1, size)]
    return gpd.GeoDataFrame({id_name: range(len(cells))}, geometry=cells, crs=buildings.crs)


def count_in(zones, id_name):
    points = buildings.assign(geometry=buildings.geometry.representative_point())
    counts = gpd.sjoin(points, zones, predicate="within").groupby(id_name).size()
    return zones[id_name].map(counts).fillna(0).astype(int)


source = grid(buildings.total_bounds, 1000, id_name="sid")
source["n"] = count_in(source, "sid")
source = source[source["n"] > 0]
target = grid(buildings.total_bounds, 700, offset=350, inset=500, id_name="tid")
truth = count_in(target, "tid")
print(f"{len(source)} source cells holding {source['n'].sum():,} buildings; "
      f"{len(target)} target cells holding {truth.sum():,}; median {int(truth.median())} per target cell")

estimates = {
    "area weighting": area_weighted(source, target, "n").values,
    "dasymetric": dasymetric(source, target, "n", buildings).values,
}
validate_areal(estimates, truth.values)

35 source cells holding 58,879 buildings; 56 target cells holding 53,137; median 722 per target cell
  source total 58,879, redistributed 50,283
  correlation between the ancillary total and n: +0.874
  source total 58,879, redistributed 58,851
  area weighting     RMSE   351.47  MAE   258.33  bias   -50.96  total     50,283 vs 53,137
  dasymetric         RMSE   458.84  MAE   247.01  bias  +102.03  total     58,851 vs 53,137
  best on RMSE: area weighting

The measurement in this article is exactly that exercise, and the sophisticated method lost.

Why area weighting is hard to beat

Why area weighting is hard to beat

It has one assumption — uniform density within each source zone — and that assumption is wrong in a way that partly cancels.

A source zone with a dense corner and an empty corner splits its value evenly between them. If the target zones average over several source zones, the over- and under-estimates partly offset, and the aggregate is close.

Measured, area weighting had a bias of −50.96 buildings per target zone against a median count of 722 — about 7%, mostly the edge loss described below. Its errors are large per zone and nearly unbiased in aggregate.

Why the dasymetric method lost here

Dasymetric weighting replaces "uniform within the zone" with "proportional to the ancillary variable". That is an improvement only if the proportionality holds.

Building footprint area against building count is not proportional. A retail park has enormous footprints and few buildings; a terrace has small footprints and many. Weighting count by area therefore moves count towards the retail park, which is wrong.

The result was a 31% worse RMSE and a positive bias of 102.03 buildings per zone.

Had the target been floor space, or population, the same ancillary variable would probably have helped substantially.

Why extensive and intensive values need different formulas

An extensive value is a total over an area, so it splits: 40% of the area carries, in expectation, 40% of the total.

An intensive value is a ratio, so it does not split. A tract with a density of 50 people per hectare contributes that density, weighted by how much of the target zone it covers — not 40% of 50.

The failure mode is silent: applying the extensive formula to a density produces a value that scales with zone size, so large target zones get implausibly large densities and the pattern is a map of zone size.

Why the totals tell you about the geometry

Area weighting redistributed 50,283 of a source total of 58,879 — a loss of 15%.

That loss is not error; it is the part of the source zones lying outside every target zone. In this measurement the target grid was offset and inset by 500 m, so the edges fell outside.

Dasymetric weighting normalised the weights within each source zone, so every source zone's full value went somewhere among the target zones it touched — preserving the total and pushing edge value inward.

Neither behaviour is correct in general. Decide which you want, and check the totals to confirm you got it.

Totals in and out, overlay cost and slivers

• Extensive values split; intensive values are area-weighted averages.
• Check the totals in and out. A large loss means partial coverage or a CRS mismatch.
• Test the ancillary variable's proportionality before using it, and validate as well — measured, it passed the correlation test and still made RMSE worse.
• Fall back to area weighting where a source zone has no ancillary features.
• Use a projected CRS. Areas in square degrees vary with latitude.
• gpd.overlay is expensive; for large layers, index first or use a spatial join on pieces.
• Small slivers from imperfect boundaries inflate the piece count; snap or filter tiny intersections.
• Report the method, the zones and the ancillary variable with the result.

Internal links

• The modifiable areal unit problem explained — why the zones decide the answer
• How to aggregate spatial data by region in GeoPandas — the simpler case of matching zones
• Overlay operations in GeoPandas: union, intersection, difference explained — the mechanics
• Spatial interpolation explained — point-based interpolation
• How to aggregate movement into flows between zones — the same zone-choice problem
• How to calculate zonal statistics in Python — raster to polygons
• How to find and fix gaps and overlaps in a polygon coverage — slivers in the overlay
• How to perform a spatial join in Python (GeoPandas) — assigning features to zones
• How to redistribute population with dasymetric mapping — using ancillary data to do better than area weights
• How to compare census years across changed boundaries — interpolation between boundary vintages

FAQ

How do I move data between two sets of polygons?

Overlay them, weight each source value by the share of the source zone in each intersection piece, and sum by target zone. That is area weighting.

What is dasymetric interpolation?

Redistributing by an ancillary variable — building footprints, land cover — rather than by area alone. It helps only when the ancillary variable is proportional to the target.

Does dasymetric weighting always beat area weighting?

No. Measured on Leiden's buildings, weighting building counts by footprint area gave an RMSE of 458.84 against area weighting's 351.47. Footprint area predicts floor space rather than count.

What is the difference between extensive and intensive values?

Extensive values (counts, totals) split between zones. Intensive values (densities, rates) do not — they are area-weighted averages.

Why do my totals not match?

Usually because source zones extend beyond the target zones, so part of their value falls outside. Check whether you want that or want the total preserved.

How do I test an ancillary variable before using it?

Correlate its total per source zone against the value you are redistributing; a correlation below about 0.7 means it may make things worse. A high one is not a guarantee — measured here, +0.874 passed and the method still lost on RMSE — so validate against known values where you can.

Should I validate areal interpolation?

Where you can — a period when both zone sets have known values. The measurement in this article is exactly that, and it overturned the expected result.