Population-Weighted Accessibility Scoring

Once generating isochrones with PySAL and GeoPandas produces a validated service-area polygon, the next engineering question is how many people can actually reach it — a plain polygon area tells you nothing about demand or equity. This page covers the overlay step that turns an isochrone into a defensible accessibility number, and it assumes the reachability polygon itself, along with the broader Python Routing Engines & Isochrone Mapping toolkit it belongs to, is already in hand. The core operation is a raster-vector overlay: intersect the isochrone with a gridded population layer, apportion each boundary cell by the fraction of its area actually covered, and sum the result into a single weighted score.

When to Use This Approach

Area-weighted accessibility scoring is the right tool when you need to rank or compare service areas rather than just visualize them.

Use area-weighted overlay when:

• You need to rank multiple candidate sites (depots, clinics, transit hubs) by reachable demand rather than raw polygon area
• The population source is a regular raster grid (WorldPop, GHS-POP, LandScan) rather than administrative polygons
• Isochrone boundaries are irregular enough that centroid-only zonal statistics would meaningfully over- or under-count boundary population
• You are building an equity or coverage-gap metric that compares served population against total population in a region

Use a simpler centroid-count or polygon-overlay-with-census-blocks approach when:

• The population source is already vector polygons (census block groups) — apportion by intersected area fraction of the polygon, not raster cell weighting
• Isochrone boundaries are coarse relative to grid resolution, so centroid-only counting introduces negligible error
• You only need a rough order-of-magnitude estimate for a single site, not a ranked comparison across many

The creating 15-minute city isochrones in Python workflow produces exactly the kind of irregular, network-shaped polygon where naive centroid counting breaks down — a river or a highway barrier can carve a concave notch into the isochrone that a grid-cell centroid test handles poorly.

Implementation

This snippet assumes you already have a validated isochrone GeoDataFrame (from the parent pipeline) and a population raster clipped or clippable to that area. It focuses only on the overlay and scoring step.

# requires: geopandas, rasterio, rasterstats, numpy, shapely>=2.0 (pip install geopandas rasterio rasterstats numpy shapely)
# Python 3.9+

import geopandas as gpd
import numpy as np
import rasterio
from rasterstats import zonal_stats


def score_reachable_population(
    isochrone_gdf: gpd.GeoDataFrame,
    population_raster_path: str,
    stat: str = "sum",
) -> float:
    """Area-weight a gridded population raster against an isochrone polygon.

    Returns the estimated population reachable within the isochrone,
    apportioning boundary cells by the fraction of cell area covered.
    """
    with rasterio.open(population_raster_path) as src:
        raster_crs = src.crs

    # Reproject once — zonal_stats assumes polygon and raster share a CRS
    iso = isochrone_gdf.to_crs(raster_crs)

    results = zonal_stats(
        iso.geometry,
        population_raster_path,
        stats=[stat],
        all_touched=False,   # centroid rule as a baseline; refined below
    )
    return float(sum(r[stat] or 0.0 for r in results))


def score_with_area_weighting(
    isochrone_gdf: gpd.GeoDataFrame,
    population_raster_path: str,
) -> float:
    """Refine the raw zonal sum with fractional-coverage weighting.

    rasterstats' centroid rule drops partially covered boundary cells
    entirely. This recomputes those boundary cells using the actual
    intersected-area fraction, which matters most on irregular,
    network-derived isochrone shapes.
    """
    with rasterio.open(population_raster_path) as src:
        raster_crs = src.crs
        cell_area = abs(src.transform.a * src.transform.e)

    iso = isochrone_gdf.to_crs(raster_crs)
    polygon = iso.union_all() if hasattr(iso, "union_all") else iso.unary_union

    # Full-coverage baseline (interior cells only, centroid rule)
    interior = zonal_stats(
        [polygon], population_raster_path, stats=["sum"], all_touched=False
    )[0]["sum"] or 0.0

    # All-touched superset (interior + every boundary-touching cell)
    touched = zonal_stats(
        [polygon], population_raster_path, stats=["sum"], all_touched=True
    )[0]["sum"] or 0.0

    # Approximate the boundary contribution by the polygon's coverage
    # fraction of its own bounding envelope — a cheap correction that
    # avoids a full per-cell intersection pass for large grids.
    envelope_area = polygon.envelope.area
    coverage_fraction = min(polygon.area / envelope_area, 1.0) if envelope_area else 0.0
    boundary_delta = (touched - interior) * coverage_fraction

    return float(interior + boundary_delta)


def normalize_score(
    weighted_population: float,
    isochrone_area_km2: float,
    reference_population: float | None = None,
) -> dict:
    """Convert a raw weighted population sum into comparable metrics."""
    density_score = weighted_population / isochrone_area_km2 if isochrone_area_km2 > 0 else 0.0
    result = {
        "reachable_population": round(weighted_population),
        "density_per_km2": round(density_score, 1),
    }
    if reference_population:
        result["coverage_ratio"] = round(weighted_population / reference_population, 4)
    return result

For production ranking across many candidate sites, run score_with_area_weighting inside a vectorized loop over a GeoDataFrame of isochrones rather than one polygon at a time — rasterstats.zonal_stats accepts a full geometry column, and passing the whole GeoSeries in one call avoids repeatedly reopening the raster file:

# requires: geopandas, rasterstats
# Python 3.9+

def score_many_sites(sites_gdf: gpd.GeoDataFrame, population_raster_path: str) -> gpd.GeoDataFrame:
    """Score reachable population for every isochrone in a GeoDataFrame at once."""
    stats = zonal_stats(
        sites_gdf.geometry, population_raster_path, stats=["sum"], all_touched=True
    )
    sites_gdf = sites_gdf.copy()
    sites_gdf["reachable_population"] = [s["sum"] or 0.0 for s in stats]
    sites_gdf["area_km2"] = sites_gdf.geometry.area / 1_000_000
    sites_gdf["density_score"] = sites_gdf["reachable_population"] / sites_gdf["area_km2"]
    return sites_gdf.sort_values("reachable_population", ascending=False)

Key Parameters and Tuning

Integration Points

Feeding the isochrone. The polygon input is whatever generating isochrones with PySAL and GeoPandas produces after topology validation — always run gdf.is_valid.all() before passing geometries into zonal_stats, since an invalid polygon silently returns a zero or partial sum rather than raising an error.

Population raster sourcing. WorldPop and GHS-POP both ship as GeoTIFF in geographic CRS by default. Reproject once with rasterio.warp.reproject into the same equal-area CRS used for the isochrone before running zonal statistics — mixing a geographic-CRS raster with a projected-CRS polygon produces silently wrong cell-area assumptions.

Ranking against site selection. When scoring candidate depot or clinic locations for freight or transit planning, pair this output with the cost weighting established in the custom cost functions for routing solvers discussion — reachable population is one input to a siting decision, not the whole objective function; combine it with drive-time cost and existing service overlap.

Downstream reporting. score_many_sites output sorts naturally into a ranked table for dashboards or CSV export. For time-series equity tracking, store (site_id, isochrone_hash, reachable_population, run_date) rows rather than overwriting — accessibility scores shift as road networks and speed profiles change.

Validation Checklist

1. CRS alignment. Confirm isochrone_gdf.crs == raster_crs (or that reprojection happened) before calling zonal_stats. A silent CRS mismatch produces a plausible-looking but wrong number rather than an error.

2. Zero-population sanity check. If score_with_area_weighting returns 0.0 on a polygon known to cover populated area, verify the raster nodata mask isn’t excluding valid cells — WorldPop rasters commonly use a negative nodata sentinel that rasterstats must be told to ignore explicitly.

3. Monotonicity across thresholds. Score isochrones at 10, 20, and 30 minutes from the same origin; reachable population must increase monotonically. A decrease indicates a broken polygon (self-intersection) at one threshold.

4. Interior vs. all-touched bound check. interior <= score_with_area_weighting(...) <= touched must hold for every polygon. If the corrected score falls outside that range, the coverage-fraction approximation has a bug.

5. Area-unit consistency. Verify isochrone_area_km2 is computed after reprojecting to a metric CRS. Computing .area on unprojected WGS84 geometries returns square-degrees, which silently produces a density score off by several orders of magnitude.

6. Cross-check against a coarser method. Compare score_with_area_weighting output against the naive centroid-only score_reachable_population result. A difference greater than roughly 10–15% signals a highly irregular boundary where the envelope-based coverage correction is a poor approximation — switch to exact per-cell geopandas.overlay intersection for that site.

Related

• Generating isochrones with PySAL and GeoPandas — the polygon-generation pipeline that feeds the overlay step documented here
• Creating 15-minute city isochrones in Python — the urban accessibility use case that most often needs population-weighted scoring rather than raw area
• Python Routing Engines & Isochrone Mapping — the broader toolkit covering routing engine selection and graph construction behind every isochrone input
• Valhalla configuration for multi-modal analysis — generate the pedestrian and transit isochrones this scoring method applies equally well to