Measuring Distance to the Nearest Protected Feature

Buffering a protected feature and testing which parcels fall inside answers one question — is this parcel within 300 feet? — and discards the number that a determination usually needs, which is how far it actually is. Buffering also scales badly: a separate buffer per threshold, recomputed whenever a threshold changes. A nearest-neighbour query returns the distance once and answers every threshold from it. This guide implements it at county scale, addresses the definitional question of what a code means by distance, and handles the parcels that legitimately find nothing. It is the distance-based sibling of the tests in proximity buffer overlap analysis.

Prerequisites

Step-by-step

Step 1: Settle what the code measures, from what, to what

The measurement is defined by two endpoints and codes differ on both. From the parcel boundary or from the proposed structure? To the feature’s edge or to its centreline? A stream setback measured to a centreline and one measured to the top of bank differ by half the channel width, which on a river is tens of feet.

ENDPOINTS = {
    "parcel_to_feature_edge":   ("parcel.geometry", "feature.geometry"),
    "structure_to_feature_edge": ("building.geometry", "feature.geometry"),
    "parcel_to_centreline":     ("parcel.geometry", "feature.centerline"),
}

Shapely’s distance returns the shortest distance between two geometries, which is boundary-to-boundary for disjoint polygons and zero for overlapping ones. Zero is a correct and important answer — a parcel containing a wetland is at distance zero from it — and code that treats zero as missing gets exactly the most constrained parcels wrong.

What the code means by distance, from what, to whatEndpoint pairs a proximity standard can name, and how far apart the answers are.Measured betweenWhere it differsParcel to feature edgeBoundary to boundaryZero when the parcel contains the featureStructure to featureedgeFootprint to boundaryNeeds a structure; often the code’s intentParcel to centrelineBoundary to a lineHalf a channel width largerAny of them,unrecordedA column called distance_ftUnreviewable — the ambiguity to remove
A stream setback measured to a centreline and one measured to the top of bank differ by half the channel width.

Step 2: Use an indexed nearest join, not a cross product

The naive implementation compares every parcel against every feature. At 300,000 parcels and 12,000 features that is 3.6 billion comparisons; the indexed version is a few seconds.

import geopandas as gpd

MAX_SEARCH_FT = 5280.0     # a mile: beyond this, "not near anything" is the answer

def nearest_feature(parcels, features, max_distance=MAX_SEARCH_FT):
    """One indexed pass: nearest feature and its distance, for every parcel."""
    joined = parcels.sjoin_nearest(
        features[["feature_id", "feature_type", "geometry"]],
        how="left",
        max_distance=max_distance,
        distance_col="distance_ft",
    )
    # sjoin_nearest emits one row per tie; keep the first deterministically.
    return (joined.sort_values(["parcel_id", "distance_ft", "feature_id"])
                  .drop_duplicates("parcel_id", keep="first"))

The tie handling is not a detail. A parcel equidistant from two features — common where features share a boundary — produces two rows, and a downstream count of “parcels within 300 feet” then over-reports. Sorting by feature identifier before deduplicating makes the choice deterministic rather than dependent on row order.

Step 3: Set the cut-off, and treat what it excludes as an answer

max_distance is what keeps the query fast, because it bounds the index search. It also creates a category of parcel with no nearest feature, and those must be recorded as “none within the search radius” rather than as null.

def resolve_distances(joined, max_distance=MAX_SEARCH_FT):
    no_hit = joined["distance_ft"].isna()
    joined.loc[no_hit, "proximity_status"] = f"none within {max_distance:.0f} ft"
    joined.loc[~no_hit, "proximity_status"] = "measured"
    joined.loc[no_hit, "distance_ft"] = float("inf")
    return joined

Using infinity rather than null is a small choice that pays off immediately: every threshold comparison downstream then works without a null check, and a parcel with no nearby feature correctly satisfies every proximity standard.

The parcels that find nothing are not missing dataA search cut-off bounds the index walk and creates a category of parcel with no nearest feature, which is an answer rather than a null.Was a protected featurefound within the searchradius?noRecord "none within the radius" and set thedistance to infinitythe parcel then satisfies every proximity standardyesRecord the distance, and the margin againsteach thresholdzero is a real answer — the parcel contains itAssert the radius exceeds the largest threshold in the packso a new threshold cannot silently exceed the search
Recording it as infinity rather than null means every threshold comparison downstream works with no null check.

Step 4: Answer every threshold from the one distance

This is the payoff over buffering. Once the distance exists, thresholds are comparisons, and adding a threshold costs nothing.

PROXIMITY_RULES = [
    {"id": "wetland-100", "feature_type": "wetland", "min_distance_ft": 100.0},
    {"id": "stream-50",   "feature_type": "stream",  "min_distance_ft": 50.0},
    {"id": "wellhead-300", "feature_type": "wellhead", "min_distance_ft": 300.0},
]

def evaluate_proximity(distances, rules=PROXIMITY_RULES):
    out = []
    for r in rules:
        subset = distances[distances["feature_type"] == r["feature_type"]]
        out.append(subset.assign(
            rule_id=r["id"],
            required_ft=r["min_distance_ft"],
            complies=subset["distance_ft"] >= r["min_distance_ft"],
            margin_ft=subset["distance_ft"] - r["min_distance_ft"]))
    return gpd.pd.concat(out, ignore_index=True)

The margin is worth carrying. A parcel at 101 feet against a 100-foot standard is compliant and worth a second look; a parcel at 900 feet is not. It is the same argument for keeping the area fraction in deciding which parcels a rule applies to — the number behind the boolean is what makes a verdict reviewable.

Step 5: Measure per feature type, not just the overall nearest

A single “nearest protected feature” column collapses distinct standards into one, and the nearest feature overall is frequently not the one that binds. Run the join per feature type.

def nearest_by_type(parcels, features, max_distance=MAX_SEARCH_FT):
    frames = []
    for ftype, group in features.groupby("feature_type"):
        near = nearest_feature(parcels, group, max_distance)
        frames.append(near.assign(feature_type=ftype))
    return gpd.pd.concat(frames, ignore_index=True)

Running one indexed join per type is also faster than it looks, because each index is smaller — and it is the arrangement that lets a single wetland layer refresh invalidate only the wetland distances.

One distance per feature type, not one overallThe nearest feature overall is usually not the one whose standard binds, so the indexed join is run per feature type and cached against that layer’s hash.Group the features by typewetland, stream, wellhead, historic siteOne indexed nearest join per typesmaller indexes, and faster than one large oneDeduplicate ties deterministicallyor every downstream count over-reportsCache on the layer hash and the endpoint ruleso a redefinition cannot reuse an old distance
Distances change when a feature layer is republished, not when the pipeline runs — so a refresh invalidates only the type that changed.

Step 6: Cache against the feature layer’s version

Distances change when a feature layer is republished, not when the pipeline runs, so recomputing them every run is waste. Key the cache on both identifiers and the layer’s content hash.

def distance_cache_key(parcel_id, feature_type, features_hash, endpoint_rule):
    return f"{parcel_id}:{feature_type}:{features_hash[:12]}:{endpoint_rule}"

Including the endpoint rule in the key is what stops a cached distance surviving a change to what the distance means — the failure that produces correct-looking numbers measured to the wrong thing.

Verification

Check that every parcel has an answer, that distances are consistent with a buffer test, and that the query really used the index.

assert distances["parcel_id"].is_unique or "feature_type" in distances, \
    "duplicate rows per parcel — the tie deduplication did not run"
assert len(distances) == len(parcels) * distances["feature_type"].nunique()

# The distance result must agree with the buffer test it replaces.
buffered = wetlands.buffer(100.0).union_all()
by_buffer = set(parcels[parcels.intersects(buffered)]["parcel_id"])
by_distance = set(distances[(distances["feature_type"] == "wetland")
                            & (distances["distance_ft"] < 100.0)]["parcel_id"])
assert by_buffer == by_distance, f"{len(by_buffer ^ by_distance)} parcels disagree"

# Zero distances are real, not missing.
inside = distances[distances["distance_ft"] == 0.0]
print(f"{len(inside)} parcels contain or touch a protected feature")

The buffer cross-check is the strongest verification available here, because it compares two genuinely different implementations of the same question. Where they disagree it is almost always the tie handling or a parcel exactly on the boundary, and both are worth seeing.

Common Pitfalls

  • A cross product instead of an indexed join. It works on a test extract and never finishes on a county.
  • Treating a zero distance as missing data. Parcels containing the feature are the most constrained ones, and they are the ones this drops.
  • Not deduplicating ties. sjoin_nearest emits a row per tied feature, which double-counts parcels in every downstream summary.
  • Measuring in a geographic CRS. Degrees are not distance, and the error varies with latitude, so the results are wrong in a way that looks orderly.
  • One distance column for all feature types. The nearest feature overall is usually not the one whose standard binds.
  • No cut-off. Without max_distance the index search is unbounded, and the parcels furthest from anything are the slowest to resolve.

Frequently Asked Questions

When is buffering still the right approach?

When the question is genuinely about area rather than distance — how much of this parcel is within the riparian buffer, which is a measurement of overlap rather than of separation. Calculating riparian buffer compliance with GeoPandas is that case. For a yes/no threshold on separation, the distance is strictly more informative and cheaper to maintain.

How does PostGIS compare here?

Favourably — this is the strongest case for the database. The <-> operator with ORDER BY and LIMIT performs an index-ordered nearest search, which is asymptotically better than a radius search that expands until it finds something. On very large feature layers the difference is substantial, and it is worth weighing against the considerations in choosing between in-database and in-process spatial joins.

Should the distance be from the parcel or from the structure?

Whichever the code says, and they are often both present — a wellhead protection standard typically measures to the proposed structure while a wetland buffer restricts the whole parcel. Compute both where both are needed and label them clearly; a column called distance_ft with no endpoint recorded is the ambiguity this whole guide is trying to remove.

How do I handle a feature layer with poor positional accuracy?

Report the distance with the layer’s stated accuracy alongside, and treat the near-threshold band as indeterminate rather than deciding it. Mapped wetland boundaries in particular are often accurate to tens of feet, which is comparable to the setbacks measured from them — so a confident verdict at 98 feet against a 100-foot standard is not supportable, and saying so is more useful than a number.

Does the cut-off bias the results?

Only by declining to measure beyond it, which is the intent. Set it comfortably above the largest threshold in the rule pack — a mile against a 300-foot standard — and assert that relationship, so that adding a larger threshold later cannot silently exceed the search radius.

Can this be tiled?

With care. Nearest-neighbour is the classic rule with unbounded reach, so a tiled run needs a halo at least as large as the cut-off, which for a mile-wide search is a large halo. In practice it is often better run as a single whole-county pass and cached, since it changes only when a feature layer changes — the exception noted in chunking county-scale runs by spatial tile.

Part of: Proximity buffer overlap analysis