Estimating Population Inside a Custom Polygon

You have a catchment polygon — a delivery zone, a fifteen-minute isochrone, a proposed store radius. Population is published for census units that do not line up with it. Areal interpolation is how you get from theirs to yours, and it is arithmetic you should do rather than buy.

Your polygon does not match theirs

Population counts are published for administrative units: census blocks and block groups in the United States, output areas in the UK, and gridded products such as Eurostat’s GEOSTAT 1 km population grid, the Global Human Settlement Layer’s GHS-POP, and WorldPop’s 100 m rasters elsewhere. Your polygon is something else entirely, so some source units fall wholly inside, some wholly outside, and the interesting ones straddle the boundary.

The whole problem is those straddling units: how much of unit i’s population to assign to your polygon. Every method is a different answer to that, and the answers differ by a lot.

Areal weighting and its one assumption

The simplest method assigns population in proportion to overlapping area:

estimate = sum over source units i of

    P_i  *  ( area(unit_i INTERSECT polygon) / area(unit_i) )

That formula assumes population is spread uniformly across each source unit. It is never exactly true and is sometimes wildly false. A census block containing a housing estate and a reservoir has all its people in part of its area; clip 40% of the block and you might have taken 0% of its residents or 100%. The size of the source unit therefore governs the size of the error: interpolating from small US census blocks is usually fine, and interpolating from a whole county is nearly meaningless.

The reassuring part is that the errors are on the boundary only. Wholly-contained units contribute their exact published count. So the uncertainty in your estimate is roughly proportional to the fraction of the total that comes from partially-overlapping units, which is a quantity you can compute and report alongside the answer.

Two implementation details bite before you get that far. The first is slivers: source-unit boundaries and your polygon are digitised from different sources, so an intersection produces hairline fragments a few square metres in size along shared edges. They contribute almost nothing to the estimate and a great deal to the runtime and the geometry-error count, so discard intersections below an area threshold. The second is invalid geometry. Self-intersecting rings and unclosed polygons are common in published boundary files, and most intersection routines either throw or return silently wrong areas on them; run a validity check and repair before the first area is computed, not after a total looks odd.

A worked interpolation

A catchment polygon overlaps three census blocks. Areas are computed in an equal-area projection, in km².

block  population   block area   intersection   share    assigned
  1        1,200        0.50           0.20        0.40       480
  2          800        0.40           0.40        1.00       800
  3        2,000        1.00           0.15        0.15       300
                                                            ------
                                        areal estimate       1,580

  fully contained contribution:                  800  (51%)
  boundary contribution:                         780  (49%)   <- the uncertainty

1,580 people, and 49% of that came from units where the uniform assumption is doing the work. That second line is the honest error bar, and it is the reason the next section exists.

Dasymetric refinement

Dasymetric mapping replaces “population is uniform over area” with “population is uniform over habitable area”, using an ancillary layer to say where people can be. Water, parks, industrial land and airports get weight zero; residential building footprints get weight proportional to their footprint area, or better, to footprint area × storeys.

Return to block 3. Its 0.15 km² intersection is mostly parkland. Suppose the block contains 40,000 m² of residential building footprint in total, and only 4,000 m² of that lies inside the polygon:

areal weight       0.15 / 1.00       = 0.15  ->  0.15 * 2,000 = 300
footprint weight   4,000 / 40,000    = 0.10  ->  0.10 * 2,000 = 200

revised total = 480 + 800 + 200 = 1,480   (down 6.3% from 1,580)

A hundred people on this polygon, and the direction of the correction is knowable in advance: clipping a park out of a block always moves the estimate down. Footprint layers are now widely available — OpenStreetMap buildings, and the open building-footprint datasets derived from imagery by footprint detection models — which makes this refinement cheap enough that there is little excuse for skipping it in built-up areas.

The projection question, answered precisely

Every step above computes areas, and area is not preserved by most projections. Compute areas in EPSG:4326 and your units are square degrees, which are not a measure of area at all. Compute them in Web Mercator, EPSG:3857, and each area is inflated by roughly 1 / cos²(latitude) — a factor of about 2.6 at 51.5°N.

Here is the part usually stated too bluntly. Areal weighting is a ratio of two areas that are metres apart, and within one census block the Mercator scale factor is essentially constant, so it cancels almost exactly in the share column. The share is not where Mercator hurts you. It hurts when you compare a share computed in Manchester with one computed in Málaga, and it hurts immediately the moment you report a density in people per km², where the constant does not cancel.

So: use an equal-area projection appropriate to the region — EPSG:5070 Albers for the conterminous US, EPSG:3035 ETRS89-LAEA for Europe, or an appropriate local equal-area or UTM zone — and do it as a matter of course rather than deciding case by case. The related trap of computing distance in the wrong reference system is a separate page, and the file-format side of getting the CRS right is in geospatial data formats.

What the number is not

1. Load the source units, reproject, and clip. Select only units whose bounding box intersects the polygon before doing exact geometry — the exact intersection is the expensive operation.
2. Compute the share per unit and sum. Report the fully-contained fraction alongside the total, as in the worked example.
3. Refine with footprints where the boundary fraction is large. If 95% of your total comes from wholly-contained units, dasymetric refinement will change the answer by almost nothing and is not worth the runtime.
4. Report residential population, and say so. Census counts are night-time, usual-residence populations. A city-centre catchment has a daytime population several times larger, and a retail-site question is almost always a daytime-population question. Confusing the two is the largest error on this page and no amount of careful geometry fixes it.