Spatial Epi Public Health GIS

Modeling Private Well Reliance in Rural Counties

This guide, part of Hazard Proximity & Drinking Water Metrics, addresses an exposure that is important, unmeasured and routinely asserted: whether a household drinks from an unregulated private well.

Problem Context & Constraints

Private wells serve a substantial rural population and fall outside drinking water regulation, so the contaminants that matter most — nitrate, arsenic, bacteria — are unmonitored for exactly the households most likely to be affected. That makes well reliance a high-value exposure variable and a difficult one, because no administrative system records it at the household level.

The naive approach is a subtraction: a household not inside a public water system service area must be on a well. It is wrong in both directions. Households outside a mapped boundary may be served by a system whose boundary was never digitised, as discussed in the sibling guide on the service area join. Households inside a boundary may still be on a well, because service areas are drawn around networks rather than around connections and rural extensions frequently pass houses that never connected.

What is available is indirect evidence. Older housing predates system expansion. Parcels distant from any distribution main are unlikely to be connected. State well-drilling registries record permits, though not abandonment. And periodic surveys have historically asked about water source at a coarse geography, giving a county or subcounty rate to calibrate against.

The constraint that makes this defensible is to build an explicit probability of well reliance per household, calibrated to whatever aggregate rate exists, and to carry that probability into the analysis rather than dichotomising it into a false certainty.

Four Weak Signals Calibrated Into One Probability Four evidence sources feed a model: housing age from parcel records, distance from the nearest distribution main or service area boundary, presence of a drilling permit in the state well registry, and lot size as a proxy for rural development pattern. Each is individually weak. The model combines them into a per-household probability, which is then calibrated so that the mean probability within each county matches the county rate from the most recent survey that asked about water source. The output is a probability, not a classification, and it is carried into the analysis as such. No single signal identifies a well; together they support a probability Housing age predates system expansion Distance to a main connection cost rises Well permit on file drilled, not necessarily used Lot size rural development pattern Per-household probability of well reliance a number between zero and one, never a yes or no Calibrate the county mean to the survey rate the only externally validated number in the chain Carry the probability into the analysis; dichotomising it discards what little is known

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, geopandas 0.14.4, scikit-learn 1.5.0
  • Parcel or address records with year built and lot size
  • Service area boundaries and, where available, distribution main geometry
  • A state well registry with permit locations
  • An aggregate well-reliance rate at county or subcounty level from a survey or census product

Step-by-Step Solution

Score each household on the available evidence, convert to a probability with a logistic form, and calibrate the within-county mean to the known aggregate rate.

# Per-household private well reliance probability, calibrated to a county rate.
# Pinned: numpy==1.26.4, pandas==2.2.2, geopandas==0.14.4, scipy==1.13.1
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
from scipy.optimize import brentq

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


def evidence_score(parcels: gpd.GeoDataFrame, mains: gpd.GeoDataFrame | None,
                   service_areas: gpd.GeoDataFrame,
                   well_permits: gpd.GeoDataFrame,
                   permit_radius_m: float = 60.0) -> pd.DataFrame:
    """Linear evidence score per parcel. Coefficients are priors, not fitted."""
    if parcels.crs.is_geographic:
        raise ValueError("project to a metric CRS before measuring distances")
    p = parcels.sort_values("parcel_id").reset_index(drop=True)

    inside = gpd.sjoin(p[["parcel_id", "geometry"]],
                       service_areas[["geometry"]], how="left", predicate="within")
    in_area = inside.groupby("parcel_id")["index_right"].first().notna()

    ref = mains if mains is not None and len(mains) else service_areas.boundary.to_frame("geometry")
    near = gpd.sjoin_nearest(p[["parcel_id", "geometry"]],
                             gpd.GeoDataFrame(geometry=ref.geometry, crs=p.crs),
                             how="left", distance_col="d_main")
    d_main = near.groupby("parcel_id")["d_main"].min()

    permit = gpd.sjoin_nearest(p[["parcel_id", "geometry"]],
                               well_permits[["geometry"]], how="left",
                               distance_col="d_permit")
    has_permit = permit.groupby("parcel_id")["d_permit"].min() <= permit_radius_m

    df = pd.DataFrame({"parcel_id": p["parcel_id"].to_numpy(),
                       "county": p["county"].to_numpy()})
    df["in_service_area"] = in_area.reindex(df["parcel_id"]).fillna(False).to_numpy()
    df["d_main_m"] = d_main.reindex(df["parcel_id"]).to_numpy()
    df["has_permit"] = has_permit.reindex(df["parcel_id"]).fillna(False).to_numpy()
    df["year_built"] = p["year_built"].to_numpy()
    df["lot_ha"] = (p.geometry.area / 10_000.0).to_numpy()

    # Prior coefficients: direction is defensible, magnitude is a judgement that
    # calibration will partly absorb. They are configuration, not fitted values.
    df["score"] = (
        1.2 * (~df["in_service_area"]).astype(float)
        + 0.9 * np.log1p(df["d_main_m"].fillna(df["d_main_m"].median()) / 500.0)
        + 1.6 * df["has_permit"].astype(float)
        + 0.5 * np.clip((1990 - df["year_built"].fillna(1990)) / 30.0, 0, 2)
        + 0.4 * np.clip(np.log1p(df["lot_ha"].fillna(0.2)), 0, 3)
    )
    log.info("evidence scored for %d parcels; %d outside any service area, "
             "%d with a well permit within %.0f m",
             len(df), int((~df["in_service_area"]).sum()),
             int(df["has_permit"].sum()), permit_radius_m)
    return df


def calibrate_to_county_rate(scored: pd.DataFrame,
                             county_rates: pd.Series) -> pd.DataFrame:
    """Shift each county's intercept so the mean probability matches the survey rate."""
    out = []
    for county, g in scored.groupby("county"):
        target = float(county_rates.get(county, np.nan))
        s = g["score"].to_numpy(float)
        if not np.isfinite(target):
            log.warning("county %s has no survey rate — probabilities left uncalibrated",
                        county)
            p = 1.0 / (1.0 + np.exp(-(s - s.mean())))
            out.append(g.assign(p_well=p, calibrated=False))
            continue

        def gap(b):
            return float(np.mean(1.0 / (1.0 + np.exp(-(s + b))))) - target

        try:
            b = brentq(gap, -30.0, 30.0)
        except ValueError:                 # target unreachable from this score spread
            b = -s.mean() + np.log(target / (1 - target))
            log.warning("county %s: exact calibration not reachable, using a "
                        "mean-matching offset", county)
        p = 1.0 / (1.0 + np.exp(-(s + b)))
        out.append(g.assign(p_well=p, calibrated=True))
        log.info("county %-14s target=%.3f achieved=%.3f offset=%+.3f n=%d",
                 county, target, float(p.mean()), b, len(g))
    return pd.concat(out, ignore_index=True)

Calibration to a county rate is what anchors the model to something externally validated. The evidence coefficients are priors — their direction is defensible and their magnitude is a judgement — so without an anchor the resulting probabilities would be an elaborate restatement of those judgements. With the anchor, the model’s job is reduced to distributing a known county total across households, which is a far weaker and more defensible claim.

Counties with no survey rate are flagged rather than given a neighbouring county’s rate. Borrowing is tempting and it silently transfers a rural county’s rate to a suburbanising one, which is the exact contrast an analysis would be trying to measure.

Calibration Distributes a Known Total, It Does Not Invent One A distribution of per-household well reliance probabilities within one county after calibration. Most households sit near zero because they are inside a service area with municipal-era housing, a smaller group sits near one because they are outside any boundary with a drilling permit on file, and a broad middle group carries intermediate probabilities. A vertical line marks the mean, which equals the county survey rate by construction. The caption records that the model decides who within the county is most likely on a well, while the survey decides how many. The survey sets how many; the model only decides which households mean = county survey rate 0 0.5 1 modelled probability of private well reliance Households in the middle carry real uncertainty and should not be forced to a side

Validation & Edge Cases

Permits without abandonment records. A drilling permit from 1974 does not mean the well is in use; many were abandoned when a system extended. Treat the permit as strong but not conclusive evidence, and where a decommissioning registry exists, use it.

Service area boundaries that overstate connection. A boundary drawn around a network includes properties along the route that never connected. Distance to a main is a better predictor than membership in a boundary where main geometry exists.

Counties with no survey rate. Leave them uncalibrated and flag them. An analysis spanning calibrated and uncalibrated counties is comparing two different quantities.

Seasonal and secondary use. Some households are connected to a system and still use a well for irrigation, or use a well in summer at a second property. The model estimates primary drinking-water source, which is what a contaminant exposure needs, but the evidence signals — a permit on file, a large lot — respond to any well at all. Where a registry records intended use, filter on it; where it does not, expect the probability to run high in areas with heavy agricultural or irrigation use and say so.

Parcel records that lag construction. Year built and lot geometry come from assessment rolls updated on their own cycle, so a subdivision completed last year may be absent or carry the parent parcel’s attributes. Compare the parcel count against a recent address point file and report the shortfall by county, because the missing parcels are concentrated exactly in the growing rural fringe where reliance is changing fastest.

Using the probability as a classification. Thresholding at 0.5 throws away the model’s entire information about uncertainty and makes the exposure look known. Where a binary is unavoidable, run the analysis across a range of thresholds and report the spread.

Newly built housing outside a boundary. Recent rural subdivisions often have shared community systems that appear in neither the service area layer nor the well registry. Flag parcels with recent construction outside all layers for manual review.

INFO water.wells evidence scored for 214880 parcels; 96412 outside any service area, 51338 with a well permit within 60 m
INFO water.wells county Marshall       target=0.412 achieved=0.412 offset=-1.884 n=18204
WARNING water.wells county Dell has no survey rate — probabilities left uncalibrated

The threshold sensitivity recommended in that list shows directly how much a dichotomised well variable depends on an arbitrary cut:

A Threshold Turns a Probability Into an Assumption Two curves plotted against the probability threshold used to classify a household as well-supplied. The share of households classified as well users falls steadily as the threshold rises from thirty to seventy per cent. The estimated association between well use and the outcome changes along with it, strengthening at higher thresholds because only the most certain households remain. The spread of estimates across plausible thresholds is the uncertainty that dichotomising hides and that carrying the probability preserves. If you must dichotomise, show what the cut point does share classified as well users estimated association 0.3 0.5 0.7 probability threshold for classifying a household as well-supplied

Compliance Notes

Record the parcel source and vintage, the service area and main geometry sources, the well registry with its coverage period, the evidence coefficients as configuration with the reasoning for each direction, the calibration target source and geography, the achieved versus target rate per county, and the list of uncalibrated counties.

Describe the variable accurately in any output: it is a modelled probability of private well reliance calibrated to an aggregate rate, not an observed water source. A study that presents it as observed overstates its precision, and one that carries the probability into the model — through weighting, or by fitting across multiple realisations drawn from it — reports an interval that reflects what is genuinely uncertain. That is usually the difference between a finding a state programme can act on and one it must first re-derive.