Spatial Epi Public Health GIS

Computing Distance to the Nearest Hazard Site at Scale

This guide, part of Hazard Proximity & Drinking Water Metrics, is about the engineering of a deceptively simple query — how far is each residence from the nearest hazardous site — when there are two million residences, forty thousand sites and several site classes.

Problem Context & Constraints

The naive implementation computes a full distance matrix. Two million residences against forty thousand sites is eighty billion distances, which is neither necessary nor feasible. A spatial index reduces it to a logarithmic search per point, and the whole job to a few minutes.

Three correctness constraints shape the implementation more than the performance one.

Distance must be metric and planar, or genuinely geodesic. Computing Euclidean distance on longitude and latitude gives a number in degrees that varies in meaning with latitude, and it is wrong by tens of per cent across a large study area. Either project to an equal-distance or local metric system and use planar distance, or use a geodesic library function. Mixing the two — projecting and then calling a geodesic function, or vice versa — produces silent nonsense.

Distance to a polygon is not distance to its centroid. A landfill is an area. The distance that matters is to its nearest edge, and for a large site the centroid can be a kilometre further away. Where the inventory gives polygons, use them; where it gives points for large facilities, note the resulting overstatement.

Per-class distances are usually what the study needs. “Nearest hazardous site” mixes a solvent recycler with a municipal landfill. Computing the nearest distance separately for each site class costs almost nothing extra with an index per class and produces a far more usable table.

Edge Distance and Centroid Distance Differ by the Site's Own Size A large irregular facility footprint with a residence to one side. The distance from the residence to the nearest point on the facility boundary is short. The distance from the residence to the facility's centroid is much longer, by roughly half the facility's width. For a small site the two are nearly identical, but for a landfill or an industrial complex spanning a kilometre the difference exceeds most buffer radii used in health studies, so the choice between them changes who is classified as exposed. For a large site, the centroid is not where the hazard is facility footprint centroid residence edge: 380 m centroid: 1 240 m With a 500 m buffer, edge distance classifies this residence exposed and centroid distance does not State which was used; the two answer different questions about a large site

Prerequisites

  • python 3.11
  • geopandas 0.14.4, shapely 2.0.4, numpy 1.26.4, pandas 2.2.2, scipy 1.13.1
  • A hazard inventory with a site class attribute and geometry — polygons where available
  • Residence points with a study identifier
  • Both layers in one projected metric coordinate reference system covering the study extent
  • The inventory extended beyond the study boundary by more than the largest distance of interest

Step-by-Step Solution

Build one index per site class, query nearest edge distance for all residences in a vectorised call, and screen the result for coordinates that behave like mailing addresses.

# Nearest-distance to hazard sites by class, for millions of residences.
# Pinned: geopandas==0.14.4, shapely==2.0.4, numpy==1.26.4, pandas==2.2.2
import logging

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

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

SUSPECT_SHARED_COORD = 25       # sites sharing a coordinate above this are suspect


def _require_planar(*frames: gpd.GeoDataFrame) -> None:
    crs = {f.crs for f in frames}
    if len(crs) != 1 or next(iter(crs)) is None:
        raise ValueError("layers must share exactly one CRS")
    if next(iter(crs)).is_geographic:
        raise ValueError("project to a metric CRS: planar distance in degrees is "
                         "not a distance")


def nearest_by_class(residences: gpd.GeoDataFrame, sites: gpd.GeoDataFrame,
                     class_col: str = "site_class") -> pd.DataFrame:
    """Nearest EDGE distance to each site class, one index per class."""
    _require_planar(residences, sites)
    residences = residences.sort_values("study_id").reset_index(drop=True)
    out = pd.DataFrame({"study_id": residences["study_id"].to_numpy()})

    for cls, grp in sites.groupby(class_col):
        grp = grp.reset_index(drop=True)
        # sjoin_nearest measures to the GEOMETRY, so polygons give edge distance.
        near = gpd.sjoin_nearest(residences[["study_id", "geometry"]],
                                 grp[["site_id", "geometry"]],
                                 how="left", distance_col="d")
        # A residence equidistant from two sites yields two rows; take the
        # smallest distance and break ties on site_id for reproducibility.
        near = (near.sort_values(["study_id", "d", "site_id"])
                    .groupby("study_id", as_index=False).first())
        out[f"dist_{cls}_m"] = near["d"].to_numpy()
        out[f"nearest_{cls}_id"] = near["site_id"].to_numpy()
        log.info("class %-18s sites=%5d median nearest=%8.0f m",
                 cls, len(grp), float(np.nanmedian(near["d"])))

    dist_cols = [c for c in out.columns if c.startswith("dist_")]
    out["dist_any_m"] = out[dist_cols].min(axis=1)
    log.info("nearest-distance complete: %d residences, %d classes",
             len(out), len(dist_cols))
    return out


def screen_site_coordinates(sites: gpd.GeoDataFrame,
                            residences: gpd.GeoDataFrame,
                            class_col: str = "site_class") -> pd.DataFrame:
    """Flag site coordinates that behave like administrative addresses."""
    _require_planar(sites, residences)
    xy = np.c_[sites.geometry.centroid.x.round(1), sites.geometry.centroid.y.round(1)]
    key = pd.Series([f"{a}_{b}" for a, b in xy])
    shared = key.value_counts()
    suspect_keys = set(shared[shared >= SUSPECT_SHARED_COORD].index)

    flags = pd.DataFrame({"site_id": sites["site_id"].to_numpy(),
                          class_col: sites[class_col].to_numpy(),
                          "shared_coordinate": key.isin(suspect_keys).to_numpy(),
                          "is_point": sites.geometry.geom_type.eq("Point").to_numpy()})
    # A "site" that many residences sit almost on top of is usually a centroid
    # fallback from the geocoder, not a facility.
    n_flag = int(flags["shared_coordinate"].sum())
    if n_flag:
        log.warning("%d site records share a coordinate with %d or more others — "
                    "likely administrative geocodes, review before use",
                    n_flag, SUSPECT_SHARED_COORD)
    return flags

Using sjoin_nearest against polygon geometry rather than centroids is what makes the distance an edge distance. It is the single line that separates a defensible metric from one that misclassifies every residence near a large facility, and it costs nothing.

The coordinate screen exists because inventories are full of administrative geocodes. When forty sites share a coordinate to within a tenth of a metre, that coordinate is a county centroid or a state fallback, not forty facilities. Those records inflate the site count near the fallback point and deflate distances there, and the pattern is easy to detect and invisible if you do not look.

Why the Index Is Not Optional Runtime plotted against the number of residences for two approaches. The full distance matrix grows linearly in the product of residences and sites, so it rises steeply and leaves the plot well before a million residences. The indexed nearest query grows close to linearly with a small constant, staying near the bottom of the plot across the whole range. A marker at two million residences shows the indexed query finishing in minutes while the matrix approach is projected in days, and a note records that memory, not time, is what actually stops the matrix approach first. Memory stops the matrix approach before the clock does full distance matrix exhausts memory here indexed nearest query 2 million residences, minutes 0 residences runtime

Validation & Edge Cases

Ties. A residence exactly between two sites yields two rows from a nearest join. Sorting on distance then site identifier makes the choice deterministic; leaving it to the join order does not.

Sites outside the study extent. If the inventory is clipped to the study boundary, every residence near the edge gets an inflated nearest distance. Extend the inventory, and verify by checking that the maximum nearest distance is not suspiciously close to a boundary-driven value.

Multipart site geometries. A facility with disjoint parcels is one record with several polygons. Edge distance handles it correctly; a centroid would fall between the parts, possibly outside the site entirely.

Zero distances. A residence inside a site polygon returns zero. That is sometimes real — housing on a former industrial parcel — and sometimes a boundary error. Count them and inspect a sample.

Class labels that are not stable. Inventory class vocabularies change between vintages. Map them to a study-defined set explicitly rather than grouping on the raw field, and record the mapping.

Distance used without direction. For an airborne hazard, undirected distance attenuates the estimate, as the parent guide describes. Where wind data exist, compute a directional variant alongside and report both.

INFO hazard.distance class landfill           sites= 1842 median nearest=   8140 m
INFO hazard.distance class solvent_recycler   sites=  614 median nearest=  19330 m
WARNING hazard.distance 312 site records share a coordinate with 25 or more others — likely administrative geocodes, review before use
INFO hazard.distance nearest-distance complete: 2118400 residences, 6 classes

The class mapping in that list is what makes per-class distances comparable across inventory vintages, and it deserves to be an explicit, versioned table:

Map Raw Classes to a Stable Study Vocabulary Raw facility class labels from two inventory vintages on the left, with differing wording and granularity: waste disposal, municipal solid waste landfill, landfill closed, solvent reclamation, and chemical recycling. Mapping lines connect them to three stable study classes on the right: landfill, solvent recycler, and other. Without the mapping, a per-class distance computed on the older vintage and one computed on the newer vintage would silently refer to different sets of facilities. Group on a declared mapping, never on the raw field raw labels, two vintages study classes waste disposal (2016) municipal solid waste landfill (2024) landfill, closed (2024) solvent reclamation (2016) chemical recycling (2024) landfill solvent recycler version the mapping table and record it with every distance run

Compliance Notes

Record the inventory identifier and extract date, the class mapping applied, whether distances were measured to polygon edges or point coordinates and for which classes, the coordinate reference system, the study-area padding, the tie-breaking rule, and the results of the coordinate screen including how many sites were flagged and what was done about them.

The computation itself consumes residence coordinates and therefore runs inside the protected environment, with an output keyed by study identifier. The site inventory and the screening results are not protected and should be published alongside any analysis: a resident or a journalist checking whether a particular facility was included in a study needs the list, and the flagged-coordinate report is exactly the kind of data-quality caveat that is better disclosed than discovered.