⚙️ How-To Interpolation & surfaces Joins & relationships Interpolation & Surfaces · step 13 of 13 Updated

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 59,391 real buildings, redistributing counts from a 1 km grid to an offset 700 m grid:

                  RMSE     MAE     bias    total
area weighting  176.69  109.82    -2.93   55,235
dasymetric      179.73  107.37   +21.66   59,391
truth                                     55,731

The dasymetric version was 1.7% worse on RMSE. That is not a coding error; it is what happens when the ancillary variable is not proportional to the target.

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.

Source zones overlaid on target zones, with each source value split by the area fraction of each intersection piece.
Overlay, weight, aggregate. The weight is where all the modelling lives.

Step-by-step solution

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 only −2.93 buildings per zone against a median zone count of 216.

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 result slightly worse — RMSE 179.73 against 176.69.

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 55,235 against a source total of 59,391
dasymetric      redistributed 59,391

Area weighting lost 4,156 because some source zones extend beyond the target grid, so part of their value falls outside every target zone. Dasymetric weighting normalised within each source zone, so it preserved the total by pushing that value into the target zones that did exist.

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.

Building footprint area weighting redistributing count badly because one large warehouse outweighs twenty houses.
Dasymetric weighting is only as good as the proportionality between the ancillary variable and the target.

Code examples

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)

The correlation test is what this article exists to argue for. Running it on the measured case would have flagged the problem before the redistribution: footprint area and building count are related but not proportional.

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.

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

Explanation

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 −2.93 buildings per target zone against a median count of 216 — under 1.5%. 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 1.7% worse RMSE and a positive bias of 21.66 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 55,235 of a source total of 59,391 — a loss of 7%.

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 slightly smaller, 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.

Four tests for an ancillary variable before using it as a dasymetric weight, and the untested shortcut that fails.
Dasymetric weighting is a hypothesis. Scoring it against area weighting is how you test it.

Edge cases or notes

  • 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 — measured, it made things 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.

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 real data, weighting building counts by footprint area gave an RMSE of 179.73 against area weighting's 176.69, because 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.

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.