How to count buildings and people exposed to a hazard zone
On this page
Problem statement
"How many people are in the flood zone?" has no single answer, because population is published for areas and hazards are mapped as their own shapes, and the two rarely nest. Every method of splitting a partly affected census unit is defensible, and they give different numbers — on the 2024 Park Fire perimeter, between 1,640 and 2,761 residents from the same census blocks, and 1,469 once population is weighted onto built-up land.
The answer is not to pick one silently. It is to compute the plausible range, use the method that best matches where people actually live, and publish the rule alongside the number.
Quick answer
Compute three apportionments — touching, centroid-inside and area-weighted — and report them together:
import geopandas as gpd, numpy as np
blocks = gpd.read_parquet("blocks_ca_north.parquet").to_crs(3310) # California Albers: equal area
fire = gpd.read_file("park_perimeter.geojson")
fire = fire[fire.poly_GISAcres > 1000].to_crs(blocks.crs)
zone = fire.union_all()
touching = blocks[blocks.intersects(zone)].copy()
touching["inside_share"] = touching.intersection(zone).area / touching.area
centroid_in = blocks[blocks.representative_point().within(zone)]
print(f"{len(touching):,} blocks touch the zone, {len(centroid_in):,} have their centre inside")
print(f"population: touching {touching.POP20.sum():,}, "
f"centroid {centroid_in.POP20.sum():,}, "
f"area-weighted {(touching.POP20 * touching.inside_share).sum():,.0f}")
print(f"housing: touching {touching.HOUSING20.sum():,}, "
f"centroid {centroid_in.HOUSING20.sum():,}, "
f"area-weighted {(touching.HOUSING20 * touching.inside_share).sum():,.0f}")320 blocks touch the zone, 224 have their centre inside
population: touching 2,761, centroid 1,677, area-weighted 1,640
housing: touching 1,381, centroid 847, area-weighted 846
Touching is an upper bound, centroid is a coin flip per block, and area-weighted assumes people are spread evenly inside a block — which in rural blocks they are not.
Fine exposure, equal-area CRS, intersection shares, weights
1. Get exposure at the finest geography available
Census blocks in the US, output areas in the UK, 100 m grids from GHSL or WorldPop elsewhere. The finer the unit, the less the apportionment rule matters.
2. Project both layers to an equal-area CRS
Area-weighting compares areas, so the CRS has to preserve them.
3. Compute the intersection share
blocks.intersection(zone).area / blocks.area is the fraction of each block inside the hazard.
4. Weight by where people are, if you can
Built-up area, building footprints or night lights redistribute population inside a block far better than area does. That is dasymetric mapping.
5. Count structures separately from people
Housing units, buildings and critical facilities each answer different questions, and a building count is not proportional to population in mixed-use areas.
6. Report the range and the rule
"1,640–2,761 residents, area-weighted estimate 1,640" says more than any single figure, and it is honest about what the data supports.
Exposure to the Park Fire with dasymetric weights
The data used on this page
Three real inputs for the 2024 Park Fire in northern California. The perimeter from the WFIGS interagency fire perimeters service (park_perimeter.geojson). The 2020 Census blocks of California with their POP20 and HOUSING20 counts, from the TIGER/Line block file (382 MB), kept for the bounding box of the main perimeter plus a 0.25° margin (blocks_ca_north.parquet) — only blocks touching the perimeter enter any result, so the margin does not change a number. And the burn severity from How to map burn severity with dNBR: its two Sentinel-2 scenes (2 July and 10 September 2024) on its 20 m window, saved as the uncorrected dNBR (dnbr.npy) and the pre-fire NBR (pre_nbr.npy).
import os
import geopandas as gpd
import numpy as np
import pystac_client
import rasterio
import requests
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
from rasterio.warp import transform_bounds
from shapely.geometry import box, shape
WFIGS = ("https://services3.arcgis.com/T4QMspbfLg3qTGWY/arcgis/rest/services/"
"WFIGS_Interagency_Perimeters/FeatureServer/0/query")
if not os.path.exists("park_perimeter.geojson"):
reply = requests.get(WFIGS, params={"where": "poly_IncidentName='PARK' AND attr_POOState='US-CA'",
"outFields": "*", "f": "geojson", "outSR": 4326}, timeout=300)
reply.raise_for_status()
open("park_perimeter.geojson", "w").write(reply.text)
if not os.path.exists("blocks_ca_north.parquet"):
fire = gpd.read_file("park_perimeter.geojson")
x0, y0, x1, y1 = fire[fire.poly_GISAcres > 1000].total_bounds # WFIGS has two small PARK fires too
blocks = gpd.read_file("https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/tl_2020_06_tabblock20.zip",
bbox=(x0 - 0.25, y0 - 0.25, x1 + 0.25, y1 + 0.25))
blocks[["GEOID20", "POP20", "HOUSING20", "ALAND20", "geometry"]].to_parquet("blocks_ca_north.parquet")
if not os.path.exists("dnbr.npy"):
BBOX = [-122.00, 40.00, -121.74, 40.22]
left, bottom, right, top = transform_bounds(4326, "EPSG:32610", *BBOX)
GRID = {"crs": "EPSG:32610", "transform": from_origin(left, top, 20, 20),
"width": int((right - left) // 20), "height": int((top - bottom) // 20)}
es = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
def scene(dates):
items = [i for i in es.search(collections=["sentinel-2-c1-l2a"], bbox=BBOX, datetime=dates,
query={"eo:cloud_cover": {"lt": 5}}).item_collection()
if shape(i.geometry).contains(box(*BBOX))]
return sorted(items, key=lambda i: i.datetime)[0]
def band(item, name):
with rasterio.open(item.assets[name].href) as src, WarpedVRT(
src, crs=GRID["crs"], transform=GRID["transform"], width=GRID["width"],
height=GRID["height"], resampling=Resampling.bilinear) as vrt:
a = vrt.read(1).astype("float32")
return np.where(a > 0, a * 0.0001, np.nan)
def nbr(item):
return (band(item, "nir08") - band(item, "swir22")) / (band(item, "nir08") + band(item, "swir22"))
pre, post = scene("2024-07-01/2024-07-10"), scene("2024-09-08/2024-09-15")
np.save("pre_nbr.npy", nbr(pre))
np.save("dnbr.npy", (nbr(pre) - nbr(post)) * 1000)
print(len(gpd.read_parquet("blocks_ca_north.parquet")), "blocks;", np.load("dnbr.npy").shape)8735 blocks; (1235, 1123)
Example 1 — where the difference comes from
partial = touching[(touching.inside_share > 0.05) & (touching.inside_share < 0.95)]
print(f"{len(partial)} blocks are partly inside; they hold {partial.POP20.sum():,} residents")
print(f"blocks fully inside hold {touching[touching.inside_share >= 0.95].POP20.sum():,}")
print("share of population in blocks that are partly inside: "
f"{partial.POP20.sum() / touching.POP20.sum():.0%}")
print(f"median touching block area: {touching.area.median() / 1e6:.2f} km², "
f"largest touching block: {touching.area.max() / 1e6:.0f} km²")74 blocks are partly inside; they hold 510 residents
blocks fully inside hold 1,271
share of population in blocks that are partly inside: 18%
median touching block area: 2.34 km², largest touching block: 111 km²
Example 2 — a dasymetric weight from land cover
import pystac_client, planetary_computer, rasterio
from rasterio.enums import Resampling
from rasterio.features import geometry_mask
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
bounds = touching.total_bounds
res = 100
grid = {"crs": touching.crs, "transform": from_origin(bounds[0], bounds[3], res, res),
"width": int((bounds[2] - bounds[0]) // res), "height": int((bounds[3] - bounds[1]) // res)}
cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace)
wc_item = next(cat.search(collections=["esa-worldcover"],
bbox=list(touching.to_crs(4326).total_bounds)).items())
with rasterio.open(wc_item.assets["map"].href) as src, WarpedVRT(
src, crs=grid["crs"], transform=grid["transform"], width=grid["width"],
height=grid["height"], resampling=Resampling.nearest) as vrt:
wc = vrt.read(1)
built = wc == 50
zone_mask = geometry_mask([zone], out_shape=(grid["height"], grid["width"]),
transform=grid["transform"], invert=True)
rows = []
for _, b in touching.iterrows():
m = geometry_mask([b.geometry], out_shape=(grid["height"], grid["width"]),
transform=grid["transform"], invert=True)
total_built = (built & m).sum()
if total_built:
rows.append(b.POP20 * (built & m & zone_mask).sum() / total_built)
else:
rows.append(b.POP20 * b.inside_share)
print(f"dasymetric estimate (population on built-up land): {sum(rows):,.0f} residents")dasymetric estimate (population on built-up land): 1,469 residents
Weighting by built-up land puts the estimate at 1,469 — below the area-weighted 1,640, because the built-up pixels in these blocks sit mostly outside the perimeter. That is the correction doing its job.
Example 3 — exposure by hazard intensity, not just inside or outside
dnbr = np.load("dnbr.npy") # Park Fire severity, 20 m
from rasterio.transform import from_origin as fo
from rasterio.warp import transform_bounds
bbox = [-122.00, 40.00, -121.74, 40.22]
l, bt, r, t = transform_bounds(4326, "EPSG:32610", *bbox)
tr = fo(l, t, 20, 20)
from shapely.geometry import box as _box
window = gpd.GeoSeries([_box(*bbox)], crs=4326).to_crs("EPSG:32610").iloc[0]
sev_blocks = blocks.to_crs("EPSG:32610")
sev_blocks = sev_blocks[sev_blocks.intersects(window)]
for lo, hi, name in [(99, 269, "low"), (269, 439, "moderate-low"), (439, 2000, "moderate-high+")]:
m = (dnbr >= lo) & (dnbr < hi)
people = 0
for _, b in sev_blocks.iterrows():
bm = geometry_mask([b.geometry], out_shape=dnbr.shape, transform=tr, invert=True)
if bm.sum():
people += b.POP20 * (m & bm).sum() / bm.sum()
print(f"{name:16} {people:6.0f} residents in the severity window")low 20 residents in the severity window
moderate-low 2 residents in the severity window
moderate-high+ 0 residents in the severity window
Twenty residents, in a window that contains 555 km² of burned forest. Most of this fire burned where nobody lived, and the exposure is concentrated at its edges — which is exactly why a single number for "people in the fire" hides more than it tells.
Why rural ranges are so wide
Why the range is so wide in rural areas
Census blocks in forested California are large — the biggest block touching this perimeter is tens of square kilometres — and a fire perimeter cuts through them at random with respect to where the houses are. A block that is 30% inside may have all its dwellings in the burned third or none of them, and the area-weighted rule assumes 30%.
Why dasymetric weighting helps
Population is not spread evenly: it sits on built-up land. Weighting each block's population by the built-up pixels inside the hazard rather than by area moves the estimate towards where people actually live. It is still a model, but a better one, and the land-cover layer is free.
Why exposure by intensity beats exposure by footprint
A fire perimeter, a floodplain or a shaking contour is a boundary drawn around a continuum. Counting people inside it throws away everything the hazard layer knows about severity. Reporting exposure per intensity class — people in moderate-high severity, dwellings in more than a metre of water — keeps it.
Daytime population, footprints and privacy
- Daytime population differs from residential population; workplaces and schools are elsewhere.
- Building footprints (Microsoft, Overture, OSM) give a better weight than land cover where available.
- Critical facilities — hospitals, schools, care homes — should be counted individually, not through a weight.
- Census vintage: 2020 blocks against a 2024 fire means four years of building.
- Privacy: block-level counts are already rounded or suppressed in some countries.
- Do not sum across hazards; the same house exposed to two hazards is one house.
Internal links
- Hazard, exposure and vulnerability explained — where exposure sits
- Exposure counts are double-counted or missed — the fix guide
- How to redistribute population with dasymetric mapping — the weighting in detail
- How to work with fire perimeter data in Python — the hazard zone used here
- How to count the population within reach of each facility — the same arithmetic for access
- Census geographies explained — blocks, tracts and why they nest
FAQ
How do I count people in a hazard zone?
Intersect the zone with census units and apportion: by touching, by centroid, by area share, or by a dasymetric weight. Report the rule.
Which method should I use?
Dasymetric where a built-up layer exists; area-weighted otherwise. Touching and centroid are bounds, not estimates.
Why do published figures differ from mine?
Different census vintage, different apportionment, or daytime rather than residential population.
Can I count buildings instead?
Yes, and it is often better — but a building count is not a population count.
What about people outside the zone who are affected?
Smoke, evacuation orders and service disruption reach much further; exposure to the footprint is a lower bound on impact.
Should I report a single number?
Report the estimate and the range from the other rules.