A point snaps to the wrong carriageway

Problem statement

A crash recorded northbound appears in the southbound records. A sign inventory has half its assets on the opposite side of the freeway. An asset that should be at postmile 58.2 comes back at 58.0. The snapping worked — it found a line, projected the point and read a measure — and it found the wrong line.

On a divided highway the two carriageways are twenty metres apart and both are plausible. Nearest-line snapping gets it right most of the time, which is the problem: the failures are rare enough to go unnoticed and frequent enough to matter.

Quick answer

Score the snap, and use the record's own side code as ground truth where it exists:

import numpy as np, pandas as pd, geopandas as gpd
from shapely.geometry import Point
from lrs import query, routes

shn = routes("County='SLO' AND Route=101")            # both carriageways
bridges = pd.DataFrame([f["attributes"] for f in
                        query("State_Highway_Bridges", "CO='SLO' AND RTE='00101'")])
bridges["PM"] = pd.to_numeric(bridges.PM, errors="coerce")
bridges = bridges.dropna(subset=["PM", "LAT", "LON"])
bridges["side"] = bridges.BRIDGE.str.strip().str[-1]  # …L and …R in the structure number
known = bridges[bridges.side.isin(["L", "R"])].copy()

pts = gpd.GeoSeries([Point(x, y) for x, y in zip(known.LON, known.LAT)],
                    crs=4326).to_crs(shn.crs)
nearest = [min(((s.geometry.distance(p), s.AlignCode) for s in shn.itertuples()),
               key=lambda t: t[0]) for p in pts]
known["offset_m"] = [n[0] for n in nearest]
known["snapped"] = [n[1][0] for n in nearest]         # 'L' or 'R'

agree = known.snapped == known.side
print(f"{len(known)} structures carry a side in their number")
print(f"nearest-line snapping agrees with it {agree.mean():.1%} of the time "
      f"({agree.sum()} of {len(known)})")
print(pd.crosstab(known.side, known.snapped).to_string())

59 structures carry a side in their number
nearest-line snapping agrees with it 96.6% of the time (57 of 59)
snapped   L   R
side           
L        28   1
R         1  29

Three per cent wrong is a small number and a large problem: the failures are not random, they cluster where the carriageways diverge, which is exactly where direction matters most.

Use the record's direction, keep and threshold the offset

1. Use the record's own direction, if it has one

A side code, a travel direction, a carriageway field. Filter the candidate routes with it before snapping, and the geometry never gets a vote.

2. Keep the offset and threshold it

The distance to the chosen line is the confidence. Snaps far from the line are the ones to look at.

3. Compare the best two candidates

If the nearest and second-nearest lines are within a few metres of each other, the choice is a coin flip. Flag it rather than recording it.

4. Snap in a projected CRS

Distances in degrees are longitude-dependent, so a "nearest" computed in 4326 is systematically wrong away from the equator.

5. Use the side of the line, not just the distance

The sign of the cross product between the line's direction and the vector to the point says which side the point lies on; combined with the route's direction it resolves most ambiguity.

6. Check the result against something independent

A side code, a direction of travel, a known postmile — anything that was not used to make the snap.

Offsets that separate right from wrong

The data used on this page

Caltrans's state highway network and its event layers — traffic counts (Traffic_AADT), bridges (State_Highway_Bridges) and functional classification — from the Caltrans ArcGIS feature services. lrs.py packages the download code from the reading guide as a module: query() pages through any layer, routes() returns centrelines with their postmiles, projected by the server to California Albers (EPSG:3310), and point_at_measure() interpolates a position along one line.

%%writefile lrs.py
"""Caltrans linear referencing: the reading guide's fetch and to_frame, as a module."""
import hashlib
import urllib.parse

import geopandas as gpd
import numpy as np
import requests
from shapely.geometry import LineString

SERVICE = "https://caltrans-gis.dot.ca.gov/arcgis/rest/services/CHhighway/{layer}/FeatureServer/0/query"
# Every page this guide requests, as Caltrans served it on 9 October 2026 (the layers are revised as
# counts are updated): FROZEN holds one file per request, named by a SHA-1 of the live URL.
FROZEN = "https://data.spatialworkflow.io/caltrans-2026-10-09/{key}.json"


def query(layer, where, out_sr=4326, want_m=False):
    """Every feature of a Caltrans layer that matches where, following the service's pages."""
    features, offset = [], 0
    while True:
        params = {"where": where, "outFields": "*", "f": "json", "outSR": out_sr,
                  "resultOffset": offset, "resultRecordCount": 2000}
        if want_m:
            params["returnM"] = "true"                   # without this there is no measure
        url = SERVICE.format(layer=layer) + "?" + urllib.parse.urlencode(params)
        key = hashlib.sha1(url.encode()).hexdigest()[:16]   # the live request: requests.get(url)
        page = requests.get(FROZEN.format(key=key), timeout=180).json()
        features += page["features"]
        if not page.get("exceededTransferLimit"):         # the service says when there is more
            return features
        offset += len(page["features"])


def to_frame(features, crs):
    """One row per path, the measures kept as an array per row."""
    rows, geoms, measures = [], [], []
    for f in features:
        for path in f["geometry"]["paths"]:
            a = np.asarray(path, dtype="float64")
            rows.append(f["attributes"])
            geoms.append(LineString(a[:, :2]))
            measures.append(a[:, 2])
    frame = gpd.GeoDataFrame(rows, geometry=geoms, crs=crs)
    frame["m"] = measures
    return frame


def routes(where, out_sr=3310):
    """State highway centrelines with postmiles, projected by the server (California Albers)."""
    frame = to_frame(query("SHN_Lines", where, out_sr=out_sr, want_m=True), out_sr)
    frame["pm_route"] = [f"{r.County}-{r.Route}{r.PMPrefix or ''}{r.PMSuffix or ''}-{r.AlignCode}"
                         for r in frame.itertuples()]         # the full route key
    return frame


def point_at_measure(line, m, target):
    """The position at a measure along one line, or None if the line does not carry it."""
    xy, m = np.asarray(line.coords), np.asarray(m)
    if not m.min() - 1e-9 <= target <= m.max() + 1e-9:
        return None
    distance = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))])
    order = np.argsort(m)                                # np.interp needs increasing x
    return line.interpolate(float(np.interp(target, m[order], distance[order])))

Writing lrs.py

Example 1 — the offset separates right from wrong

print(f"offset when the snap agrees with the side code:  "
      f"median {known.offset_m[agree].median():5.1f} m, "
      f"90th {known.offset_m[agree].quantile(0.9):5.1f} m")
print(f"offset when it disagrees:                        "
      f"median {known.offset_m[~agree].median():5.1f} m")
print(known[~agree][["BRIDGE", "NAME", "PM", "side", "snapped", "offset_m"]]
      .round(1).to_string(index=False))

for cut in (10, 20, 30):
    flagged = known.offset_m > cut
    caught = (flagged & ~agree).sum()
    print(f"flagging snaps over {cut:2} m: {flagged.sum():2} flagged, "
          f"{caught} of {(~agree).sum()} errors caught")

offset when the snap agrees with the side code:  median   5.4 m, 90th   7.9 m
offset when it disagrees:                        median  27.7 m
  BRIDGE                    NAME   PM side snapped  offset_m
49 0162L     SOUTH SAN MIGUEL UC 65.1    L       R      19.3
49 0185R LOS OSOS VALLEY ROAD OC 25.9    R       L      36.1
flagging snaps over 10 m:  4 flagged, 2 of 2 errors caught
flagging snaps over 20 m:  1 flagged, 1 of 2 errors caught
flagging snaps over 30 m:  1 flagged, 1 of 2 errors caught

A ten-metre threshold catches both errors at the cost of flagging four records — a review list of four instead of a silent error rate of two. Note that the two wrong snaps are both overcrossings, structures that carry another road over this one and whose published coordinate sits wherever the crossing road happens to be.

Example 2 — which side of the line the point is on

def side_of(line, point):
    """+1 left of the direction of travel, -1 right, using the nearest segment."""
    xy = np.asarray(line.coords)
    d = line.project(point)
    along = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))])
    i = min(np.searchsorted(along, d), len(xy) - 1)
    a, b = xy[max(i - 1, 0)], xy[min(i, len(xy) - 1)]
    cross = ((b[0] - a[0]) * (point.y - a[1]) - (b[1] - a[1]) * (point.x - a[0]))
    return 1 if cross > 0 else -1

right_only = shn[shn.AlignCode == "Right"]
sides = []
for p, s in zip(pts, known.side):
    seg = min(right_only.itertuples(), key=lambda r: r.geometry.distance(p))
    sides.append((side_of(seg.geometry, p), s))
tab = pd.crosstab(pd.Series([x[0] for x in sides], name="geometric side"),
                  pd.Series([x[1] for x in sides], name="side code"))
print(tab.to_string())
print("\nthe sign is consistent within a carriageway, which is what makes it usable "
      "as a tie-break rather than as a label")

side code        L   R
geometric side        
-1               0  29
 1              29   1

the sign is consistent within a carriageway, which is what makes it usable as a tie-break rather than as a label

Example 3 — filter the candidates instead of ranking them

def snap_with_side(point, side, frame, max_offset=60.0):
    candidates = frame[frame.AlignCode.str[0] == side]
    if not len(candidates):
        return None, np.nan
    best = min(((r.geometry.distance(point), r) for r in candidates.itertuples()),
               key=lambda t: t[0])
    return (best[1], best[0]) if best[0] <= max_offset else (None, best[0])

results = [snap_with_side(p, s, shn) for p, s in zip(pts, known.side)]
ok = [r for r in results if r[0] is not None]
print(f"{len(ok)} of {len(results)} snapped within 60 m of the correct carriageway")
print(f"agreement with the side code: 100% by construction")
offsets = np.array([r[1] for r in results if r[0] is not None])
print(f"offsets: median {np.median(offsets):.1f} m, max {offsets.max():.1f} m")
rejected = [(b.BRIDGE, round(r[1], 1)) for b, r in zip(known.itertuples(), results)
            if r[0] is None]
print(f"rejected as too far: {rejected}")

59 of 59 snapped within 60 m of the correct carriageway
agreement with the side code: 100% by construction
offsets: median 5.4 m, max 56.9 m
rejected as too far: []

Nothing is rejected here and the median offset is unchanged, but the maximum rises from 36 to 57 metres: those are the records that were previously snapping to the wrong line at a comfortable-looking distance. Filtering by side first turns an accuracy problem into a coverage problem, which is the better kind — the results that survive are right, and anything that cannot be matched is visible instead of silently mislocated.

Nearest is right often enough to be dangerous

Why nearest-line is right often enough to be dangerous

Two carriageways twenty metres apart, and assets that sit beside one of them, produce a correct nearest answer almost every time. The exceptions are structures that span both, points digitised from the median, and places where the alignments cross — under an interchange, at a bridge, through a town. Those are precisely the locations that matter for a safety or maintenance analysis.

Why the offset is the best available diagnostic

A wrong snap is usually a distant snap, because the point was really beside the other line. In the Caltrans structures here, the agreeing snaps sit a few metres from their carriageway and the disagreeing ones several times further. That ratio is what makes a threshold work — not because far snaps are always wrong, but because they are the ones worth looking at.

Why direction beats geometry

The record usually knows something the geometry does not: which way the traffic was going, which side the sign faces, which carriageway the crew worked on. Using that to filter candidates before snapping removes the ambiguity entirely, and leaves the geometry to do the one job it can do reliably — finding the position along the line that was already chosen.

Spanning structures, ramps and median assets

• Structures that span both carriageways have one record and two valid answers; the side code usually says which is primary.
• Independent alignments — where the two directions take different streets — break the twenty-metre assumption entirely.
• Ramps are closer to a point than either mainline in the middle of an interchange.
• Median assets are equidistant by design; use the direction they serve.
• Three-line cases: HOV or express lanes add a third candidate on the same corridor.
• Batch runs should record the second-best candidate and its offset, so ambiguity is auditable afterwards.

Internal links

• How to find the measure of a point on a route — the operation this fixes
• How to transfer events between two route systems — moving events between carriageways deliberately
• Postmile and milepost systems explained — where the alignment code comes from
• Linear referencing explained: routes, measures and events — the route key
• How to snap points to a street network in Python — the general version
• Street network routing returns absurd distances — the routing cousin of this bug

FAQ

Why did my point snap to the opposite carriageway?

The other line was nearer at that spot — usually near an interchange, a structure or a median asset.

How do I stop it?

Filter candidate routes by the record's direction or side code before snapping, rather than ranking by distance.

What offset should I allow?

Enough to reach the correct carriageway and not the other one: on a rural divided highway, tens of metres.

How do I detect wrong snaps without ground truth?

Compare the best and second-best candidates; a small margin means the choice was arbitrary.

Does projecting in the right CRS matter?

Yes — nearest computed in degrees is not nearest on the ground.

What if the record has no direction?

Record the ambiguity with the result, and treat direction-specific analyses of that record as unsupported.