Spatial Epi Public Health GIS

Zonal Statistics for Raster Exposure at Address Points

This guide, part of Assigning Exposure to Residential Locations, covers the mechanics of getting values out of an exposure raster at scale — the step that looks trivial, runs for six hours, and quietly returns a half-cell offset if the alignment is not checked.

Problem Context & Constraints

The obvious implementation opens the raster, loops over addresses, and reads a window per record. For a cohort of a million it performs a million seeks, and if the raster is remote it performs a million range requests. The same work done as a vectorised read of the whole array, or a blocked read of the tiles that actually contain addresses, finishes in minutes.

Performance, though, is the least dangerous part. Three correctness issues cost more.

Coordinate reference system mismatch. Reprojecting the points to the raster’s system is right; reprojecting the raster to the points’ system resamples the data and changes the values. The direction matters, and doing it the wrong way round is a common and invisible error.

Cell-centre versus cell-corner indexing. A geotransform maps a coordinate to a cell index, and an implementation that treats the transform origin as a cell centre rather than a corner shifts every value by half a cell. On a 1 km grid that is 500 m of systematic error applied to the whole cohort, in one direction.

Nodata that is not declared. A raster carrying −9999 as a fill value without declaring it in the header returns −9999 as an exposure. It will not raise; it will produce a mean that is obviously wrong only if somebody looks.

The constraint that ties these together is that the sampling code must assert what it assumes rather than trusting the file’s metadata, because exposure rasters are frequently assembled by hand from sources that disagree with one another.

The Half-Cell Offset and What It Costs A grid of raster cells with an address point placed inside one of them. In the correct indexing, the transform origin is the outer corner of the first cell and the point falls in the cell whose value is twelve. In the offset indexing, the origin is treated as the centre of the first cell, shifting the grid by half a cell so the same point now falls in a neighbouring cell whose value is nineteen. Because the shift applies to every record in the same direction, the error is systematic rather than random and does not average out across a cohort. A half-cell shift is systematic, so it never averages out Correct: origin is a cell corner 8 19 21 9 12 17 7 10 14 sampled value = 12 Offset: origin treated as a centre 8 19 21 9 12 17 7 10 14 sampled value = 19 Same coordinate, same raster, different value — and the error has a consistent sign Test the indexing against a cell whose value you know before running the cohort

Prerequisites

  • python 3.11
  • rasterio 1.3.10, numpy 1.26.4, geopandas 0.14.4, pandas 2.2.2
  • An exposure raster with a declared coordinate reference system and, ideally, a declared nodata value
  • Geocoded address points with a study identifier and geocode match type
  • Enough memory for one raster band, or a tiled access pattern where there is not

Step-by-Step Solution

Reproject points to the raster, convert coordinates to indices in one vectorised operation, and read the array once. Verify the transform before anything else.

# Vectorised raster sampling at address points, with alignment and nodata assertions.
# Pinned: rasterio==1.3.10, numpy==1.26.4, geopandas==0.14.4, pandas==2.2.2
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("exposure.zonal")


def verify_raster(path: str, known_point: tuple[float, float] | None = None,
                  known_value: float | None = None) -> dict:
    """Assert what the sampling code is about to assume about this raster."""
    with rasterio.open(path) as src:
        t = src.transform
        meta = {"crs": src.crs.to_string() if src.crs else None,
                "res_x": abs(t.a), "res_y": abs(t.e),
                "nodata": src.nodata, "dtype": src.dtypes[0],
                "width": src.width, "height": src.height}
        if meta["crs"] is None:
            raise ValueError("raster has no CRS — refusing to sample")
        if meta["nodata"] is None:
            log.warning("raster declares no nodata value; sentinel fills will be "
                        "returned as real exposures")
        if abs(t.b) > 1e-12 or abs(t.d) > 1e-12:
            raise ValueError("rotated geotransform — reproject the raster first")
        if known_point is not None:
            # Ground-truth the indexing: one cell whose value you verified by hand.
            row, col = src.index(*known_point)
            got = float(src.read(1, window=((row, row + 1), (col, col + 1)))[0, 0])
            if known_value is not None and not np.isclose(got, known_value):
                raise AssertionError(
                    f"indexing check failed: expected {known_value} at "
                    f"{known_point}, read {got} — likely a half-cell offset")
            log.info("indexing check passed at row=%d col=%d value=%.4f", row, col, got)
    log.info("raster verified: %s", meta)
    return meta


def sample_points(path: str, points: gpd.GeoDataFrame, band: int = 1) -> pd.DataFrame:
    """Sample a raster at every point in one vectorised pass."""
    points = points.sort_values("study_id").reset_index(drop=True)
    with rasterio.open(path) as src:
        pts = points.to_crs(src.crs) if points.crs != src.crs else points
        xs = pts.geometry.x.to_numpy()
        ys = pts.geometry.y.to_numpy()
        rows, cols = rasterio.transform.rowcol(src.transform, xs, ys)
        rows = np.asarray(rows)
        cols = np.asarray(cols)
        inside = (rows >= 0) & (rows < src.height) & (cols >= 0) & (cols < src.width)
        arr = src.read(band, masked=True)
        vals = np.full(len(pts), np.nan)
        vals[inside] = arr[rows[inside], cols[inside]].filled(np.nan)
        nod = src.nodata

    out = pd.DataFrame({"study_id": points["study_id"].to_numpy(), "exposure": vals})
    if nod is not None:
        out.loc[np.isclose(out["exposure"], nod), "exposure"] = np.nan
    # Sentinel sweep: values far outside a plausible range are fills, declared or not.
    suspicious = out["exposure"].abs() > 1e5
    if suspicious.any():
        log.error("%d sampled values look like undeclared fills — inspect before use",
                  int(suspicious.sum()))
        out.loc[suspicious, "exposure"] = np.nan
    log.info("sampled %d points: outside_extent=%d missing=%d",
             len(out), int((~inside).sum()), int(out["exposure"].isna().sum()))
    return out

The verify_raster call is the part worth insisting on. It takes one known coordinate with a value you checked by hand in a desktop tool, and it fails loudly if the indexing does not agree. Every other check in the function — rotation, declared coordinate system, nodata — protects against a class of file that arrives in exposure work regularly, assembled from mosaics and reprojected by somebody else.

For buffer means rather than point lookups, the equivalent vectorised approach is to precompute a cumulative sum over the array and read each buffer’s rectangular approximation in constant time. The approximation is a square rather than a circle, which is acceptable when the buffer is small relative to the surface’s spatial correlation and unacceptable when it is not — a decision to record, not to make silently.

Three Access Patterns at Cohort Scale Three sampling strategies compared for a cohort of one million address points. Per-record window reads perform one seek per record and scale badly, particularly against remote storage. A single vectorised array read converts all coordinates to indices at once and reads the band once, which is the right approach whenever the band fits in memory. An integral-image path precomputes a cumulative sum so that any rectangular buffer mean is available in constant time, at the cost of approximating circular buffers with squares. Each option lists the situation it suits and the cost it carries. Choose the access pattern before optimising anything else Per-record window one seek per address exact for any geometry suits: a few thousand records cost: hours against remote storage Vectorised array read all indices computed at once band read exactly once suits: point lookups at any scale cost: one band must fit in memory Integral image cumulative sum precomputed any rectangle in constant time suits: many buffer means cost: squares approximate circles Whichever path is used, verify the indexing against a hand-checked cell first an offset costs the same whether it is computed quickly or slowly

Validation & Edge Cases

Points outside the raster extent. These return no value and must be counted, not silently nulled. A block of them usually means a coordinate reference system problem rather than genuinely distant addresses.

Masked arrays and the fill value. Reading with masked=True and then calling filled without an explicit fill returns the array’s default fill, which is a large number. Always pass the fill explicitly.

Integer rasters carrying scaled values. Land surface temperature is often stored as scaled integers. Applying the scale factor and offset is the caller’s job in most libraries, and skipping it produces exposures in the tens of thousands.

Multi-band rasters with band order assumptions. Prediction in band 1 and variance in band 2 is a convention, not a guarantee. Read the band descriptions and assert them.

Duplicate coordinates. Apartment buildings give many records the same value, which is correct. But a sudden spike of identical coordinates that is not a residential building is a geocoder fallback, and its sampled exposure is meaningless.

INFO exposure.zonal indexing check passed at row=2184 col=917 value=11.4200
INFO exposure.zonal raster verified: {'crs': 'EPSG:5070', 'res_x': 1000.0, 'res_y': 1000.0, 'nodata': -9999.0, 'dtype': 'float32', 'width': 4881, 'height': 3109}
INFO exposure.zonal sampled 1042118 points: outside_extent=214 missing=1908

The coordinate reference system direction in that list is the error with the least obvious symptom, so it helps to see both paths:

Move the Points, Never the Raster Two paths for sampling an exposure raster in a different coordinate system from the addresses. In the correct path the address points are transformed into the raster's system, an exact operation, and each point reads an original cell value. In the wrong path the raster is warped into the points' system, which resamples every cell and blends neighbouring values, so the sampled exposure is an interpolated number that never existed in the source surface. Both paths run without error. Transforming points is exact; warping a raster changes its values Points to raster system coordinate transform only raster untouched reads original cell values correct Raster to points system every cell resampled neighbours blended reads interpolated values wrong, and silent Both run without error; only the first returns numbers that exist in the surface

Compliance Notes

Record the raster identifier, vintage, band descriptions, coordinate reference system, resolution, declared nodata and the result of the indexing verification including the coordinate and value used. Record the access pattern and, where an integral-image approximation was used for buffers, the fact that circular buffers were approximated by squares and at what radius.

The sampling step runs inside the protected environment because it consumes address coordinates. The log must carry counts and raster metadata only — never the coordinate that failed, never the study identifier of a record outside the extent. Where a verification coordinate is hard-coded for the indexing check, use a public landmark rather than a case address.