Spatial Epi Public Health GIS

Linking Case Addresses to Public Water System Boundaries

This guide, part of Hazard Proximity & Drinking Water Metrics, covers the join that sits under every drinking water epidemiology study: deciding which regulated system, if any, supplied the water at a case address.

Problem Context & Constraints

A point-in-polygon join looks like the whole task and is the easy part. The difficulties are in the boundary data.

Boundaries overlap. A wholesale system’s service area frequently contains the retail systems that purchase from it, so an address falls inside two or three polygons that are all correct. A naive join returns whichever row the spatial index yielded first, which is arbitrary and irreproducible across library versions.

Boundaries are approximate. Many published service areas are drawn from billing ZIP codes, from parcel lists of unknown vintage, or from a utility’s own sketch. A point-in-polygon test against such a polygon returns a crisp answer from a fuzzy input, and the crispness is misleading near edges.

Boundaries are missing. For large parts of many states no digitised service area exists at all. Addresses there are not unserved; they are unmapped, and treating a missing boundary as “on a private well” converts a data gap into an exposure classification.

The constraint that makes this tractable is to treat the join as producing an assignment plus a confidence, rather than a bare system identifier. Downstream analyses can then restrict by confidence, and the fraction of records at each level becomes a reportable property of the study rather than an invisible assumption.

Overlapping Service Areas and the Resolution Rule A large wholesale system service area containing two smaller retail service areas, with a gap between them that belongs to the wholesaler directly. Three addresses are marked: one inside a retail area and therefore inside two polygons, one in the gap and therefore inside only the wholesale area, and one outside every polygon. The resolution rule assigns the first address to the retail system because it is the smaller and more specific polygon, the second to the wholesale system directly, and the third to no system with a confidence flag rather than to a private well. Inside two polygons is the normal case, not an error wholesale system service area retail system A retail system B in two polygons wholesale only outside everything unmapped, not unserved assign no system and flag never infer a private well Rule: smallest containing polygon wins, ties broken on a stable system identifier

Prerequisites

  • python 3.11
  • geopandas 0.14.4, pandas 2.2.2, numpy 1.26.4
  • Service area polygons with a system identifier, a system type (community, wholesale, transient) and a boundary provenance attribute where available
  • Geocoded case addresses with a match type
  • A purchase relationship table where wholesale supply is to be traced

Step-by-Step Solution

Join, resolve overlaps by the smallest containing polygon, and emit a confidence level per record derived from both the geocode quality and the boundary provenance.

# Deterministic address-to-water-system join with a per-record confidence.
# Pinned: geopandas==0.14.4, pandas==2.2.2, numpy==1.26.4
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("water.join")

GOOD_GEOCODES = {"rooftop", "parcel"}
GOOD_BOUNDARIES = {"surveyed", "parcel_derived"}
EDGE_BUFFER_M = 150      # a match this close to a boundary is not confident


def assign_water_system(cases: gpd.GeoDataFrame, systems: gpd.GeoDataFrame,
                        geocode_col: str = "geocode_match",
                        provenance_col: str = "boundary_provenance") -> pd.DataFrame:
    """Assign each address to a system, resolving overlaps and grading confidence."""
    if cases.crs != systems.crs or cases.crs.is_geographic:
        raise ValueError("both layers must share one projected metric CRS")
    cases = cases.sort_values("study_id").reset_index(drop=True)
    sys_ = systems.copy()
    sys_["poly_area"] = sys_.geometry.area

    hits = gpd.sjoin(cases[["study_id", geocode_col, "geometry"]],
                     sys_[["pwsid", "system_type", provenance_col,
                           "poly_area", "geometry"]],
                     how="left", predicate="within")

    # Overlap resolution: the SMALLEST containing polygon is the most specific
    # supplier. Ties break on the system identifier so the result is stable.
    hits = hits.sort_values(["study_id", "poly_area", "pwsid"])
    n_multi = int((hits.groupby("study_id").size() > 1).sum())
    chosen = hits.groupby("study_id", as_index=False).first()

    # Distance to the boundary of the chosen polygon: an address 20 m inside an
    # approximate line is not confidently inside it.
    geom = cases.set_index("study_id").geometry
    edge_dist = []
    sys_by_id = sys_.set_index("pwsid").geometry
    for rec in chosen.itertuples():
        if not isinstance(rec.pwsid, str):
            edge_dist.append(np.nan)
            continue
        edge_dist.append(float(geom.loc[rec.study_id].distance(
            sys_by_id.loc[rec.pwsid].boundary)))
    chosen["dist_to_boundary_m"] = edge_dist

    def grade(r):
        if not isinstance(r.pwsid, str):
            return "unassigned"
        if (r._asdict()[geocode_col] in GOOD_GEOCODES
                and r._asdict()[provenance_col] in GOOD_BOUNDARIES
                and (r.dist_to_boundary_m or 0) >= EDGE_BUFFER_M):
            return "high"
        if r._asdict()[geocode_col] in GOOD_GEOCODES:
            return "medium"
        return "low"

    chosen["assignment_confidence"] = [grade(r) for r in chosen.itertuples()]
    counts = chosen["assignment_confidence"].value_counts().to_dict()
    log.info("water system join: %d addresses, %d fell in multiple polygons, "
             "confidence=%s", len(chosen), n_multi, counts)
    return chosen[["study_id", "pwsid", "system_type", "dist_to_boundary_m",
                   "assignment_confidence"]]


def trace_purchases(assigned: pd.DataFrame, purchases: pd.DataFrame,
                    max_depth: int = 3) -> pd.DataFrame:
    """Follow seller relationships so a household inherits its wholesaler too.

    purchases: buyer_pwsid, seller_pwsid
    """
    chain = {r.buyer_pwsid: r.seller_pwsid for r in purchases.itertuples()}
    out = []
    for rec in assigned.itertuples():
        ids, cur, depth = [rec.pwsid], rec.pwsid, 0
        while cur in chain and depth < max_depth:
            cur = chain[cur]
            if cur in ids:                    # a reported cycle; stop and flag
                log.warning("purchase cycle at %s — chain truncated", cur)
                break
            ids.append(cur)
            depth += 1
        out.append({"study_id": rec.study_id, "pwsid_chain": ids,
                    "chain_depth": len(ids) - 1})
    df = pd.DataFrame(out)
    log.info("purchase tracing: median chain depth=%.0f, %d with a wholesaler",
             df["chain_depth"].median(), int((df["chain_depth"] > 0).sum()))
    return df

Resolving overlaps by smallest containing polygon is a rule, not a heuristic: the retail system that directly bills a household is the more specific supplier, and it is the one whose distribution-system violations — lead and copper, disinfection byproducts — describe the water at the tap. The wholesaler is reached separately through the purchase chain, which is why both are kept.

The distance-to-boundary check is what converts a crisp join into an honest one. An address 20 m inside an approximately drawn boundary is not confidently inside it, and grading those to medium confidence lets a sensitivity analysis exclude them without discarding the whole dataset.

Two Inputs Decide the Assignment Confidence A grid with geocode quality on one axis and boundary provenance on the other. A rooftop or parcel geocode combined with a surveyed or parcel-derived boundary, and comfortably away from the boundary line, gives high confidence. A good geocode with an approximate boundary gives medium confidence, as does a good geocode close to a boundary line. A coarse geocode gives low confidence regardless of boundary quality, because the address could be on either side. The grid shows that the weaker of the two inputs sets the grade. The weaker input sets the grade, and the grade travels with the record high and not near a boundary line medium boundary is a sketch low a coarse geocode could be on either side of any boundary surveyed or parcel-derived boundary approximate boundary rooftop or parcel coarse geocode Report the share at each grade; restricting to high confidence is a valid sensitivity analysis

Validation & Edge Cases

Addresses in a gap between adjacent systems. Real service areas tessellate imperfectly, leaving slivers. An address in a sliver is almost certainly served by one of the neighbours; assign to the nearest system within a stated distance and grade it low, rather than declaring it unassigned.

Systems with multiple disjoint polygons. A utility may serve several separate towns under one identifier. The join handles it naturally, but area-based overlap resolution can then pick a large multipart wholesale polygon over a small retail one. Compare the area of the specific part, not the multipart total.

Purchase cycles. Reciprocal supply agreements produce loops in the purchase table. Cap the traversal depth and log the cycle rather than looping forever.

Boundary vintage against case date. Service areas expand. A 2024 boundary applied to a 2011 case can place a household inside a system it did not join until 2018. Where dated boundaries exist, use them; where not, record the vintage mismatch.

Apartment and mobile-home park addresses. A large building or park is frequently served through a single master meter, and in some states the park itself is a regulated system with its own identifier. A join against the surrounding municipal boundary then assigns the wrong supplier to every household in it. Where a master-metered system exists in the registry, its polygon is often tiny or absent, so check the registry by name as well as by geometry for known multi-unit addresses.

Systems that changed identifier. Consolidations and transfers reissue system identifiers, so a violation history keyed to an old identifier will not join to a boundary keyed to the new one. Maintain an identifier crosswalk with effective dates, and report how many records were matched through it rather than directly.

Assuming unassigned means private well. This is the most consequential error in the whole join. Unassigned means unmapped. Where private well reliance is the exposure of interest, it needs its own evidence — a survey, a state well registry, a modelled estimate — and its uncertainty reported separately.

INFO water.join water system join: 51204 addresses, 18886 fell in multiple polygons, confidence={'high': 29118, 'medium': 12904, 'low': 6141, 'unassigned': 3041}
INFO water.join purchase tracing: median chain depth=0, 14220 with a wholesaler

The sliver case in that list arises wherever two utilities digitised their shared boundary separately, and the nearest-system rule handles it without inventing a well:

An Address in the Gap Between Two Utilities Two adjacent water system service areas whose shared boundary was digitised separately, leaving a thin sliver between them that belongs to neither polygon. An address falls in the sliver. A point-in-polygon join assigns it to no system, which a naive rule would treat as a private well. The correct handling assigns it to the nearest system within a stated distance, grades the assignment low confidence, and records the distance, so the record stays in the analysis without being mistaken for an unregulated supply. In the gap is not the same as off the network system A system B address in sliver naive: no system called a well, wrongly nearest: system A low confidence, 4 m assign to the nearest system within a stated distance, grade it low, record the distance

Compliance Notes

Record the boundary dataset identifier and vintage, the provenance attribute definitions, the overlap resolution rule, the edge buffer distance, the confidence grading rules, the counts at each confidence level, the purchase table source and the traversal depth cap.

The join runs inside the protected environment because it consumes case addresses; the output is keyed by study identifier and carries a system identifier, which is not itself identifying. Report the confidence distribution in any publication, and run the primary analysis restricted to high-confidence assignments as a sensitivity check. A water-quality association that appears only when low-confidence assignments are included is describing the geocoding, not the water.