How to aggregate events to route segments

On this page

Problem statement

Point events on a route — structures, signs, collisions, complaints — are counted per stretch of road before they mean anything. The choice of stretch decides the answer: fixed miles hide anything shorter than a mile, segments from another table put boundaries in arbitrary places, and a sliding window finds concentrations that a fixed grid splits in half.

This guide counts the 93 structures on US 101 in San Luis Obispo county — from an inventory of 193 across the county — three ways, and turns one of them into a rate per vehicle-mile.

Quick answer

Bin by measure, and be explicit about the bin:

import numpy as np, pandas as pd
from lrs import query

bridges = pd.DataFrame([f["attributes"] for f in query("State_Highway_Bridges", "CO='SLO'")])
bridges["PM"] = pd.to_numeric(bridges.PM, errors="coerce")
bridges["LENG"] = pd.to_numeric(bridges.LENG, errors="coerce")
on101 = bridges[(bridges.RTE == "00101")].dropna(subset=["PM"])

edges = np.arange(0, 70.0 + 1e-9, 1.0)                 # one-mile bins
counts, _ = np.histogram(on101.PM, bins=edges)
busiest = counts.argmax()
print(f"{len(on101)} structures on US 101, postmiles "
      f"{on101.PM.min():.2f}–{on101.PM.max():.2f}")
print(f"{(counts > 0).sum()} of {len(counts)} miles carry at least one")
print(f"busiest mile: {edges[busiest]:.0f}–{edges[busiest + 1]:.0f} "
      f"with {counts[busiest]} structures")
print(f"mean {counts.mean():.2f} per mile, max {counts.max()}")

93 structures on US 101, postmiles 0.09–67.23
38 of 70 miles carry at least one
busiest mile: 16–17 with 10 structures
mean 1.33 per mile, max 10

A one-mile bin is a choice, not a neutral default: it decides whether two structures 200 feet apart on either side of postmile 12 are one cluster or two single miles.

The same events, two binnings: the fixed grid splits the cluster the window finds.

One route, a bin for the question, weights and exposure

1. Filter to one route and one measure system

Postmiles from several routes in one histogram are meaningless — they all start at zero.

2. Choose the bin for the question

Fixed bins for reporting, existing segments for joining to other attributes, a sliding window for finding concentrations.

3. Weight by whatever the event carries

A count treats a 500-foot viaduct and a culvert alike. Summing deck area, lane-miles or cost usually answers the real question.

4. Normalise by exposure, not just by length

Events per mile compares stretches of road. Events per vehicle-mile compares risk, and needs the volume from another table.

5. Watch the bins with no road in them

A bin inside an equation gap has no road, so a zero there is not a zero rate.

6. Report the bin definition with the number

"Three per mile" means nothing without the bin width, the route and the measure range.

Per mile and per vehicle-mile rank the same road differently.

Three binnings and a rate per mile

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 — three binnings of one dataset

from lrs import routes

shn = routes("County='SLO' AND Route=101 AND AlignCode='Right'").sort_values("bPM")

# a) fixed one-mile bins
fixed = pd.DataFrame({"from_pm": edges[:-1], "to_pm": edges[1:], "n": counts})

# b) the route's own segments
seg_rows = []
for s in shn.itertuples():
    inside = on101[(on101.PM >= s.bPM) & (on101.PM < s.ePM)]
    seg_rows.append({"from_pm": s.bPM, "to_pm": s.ePM, "miles": s.ePM - s.bPM,
                     "n": len(inside), "deck_ft": inside.LENG.sum()})
by_segment = pd.DataFrame(seg_rows)

# c) a sliding window
step, width = 0.1, 1.0
starts = np.arange(0, 69.0 + 1e-9, step)
window = np.array([((on101.PM >= a) & (on101.PM < a + width)).sum() for a in starts])
peak = window.argmax()

print(f"fixed miles      : {len(fixed):3} bins, max {fixed.n.max()} structures")
print(f"route segments   : {len(by_segment):3} bins, max {by_segment.n.max()} structures "
      f"in {by_segment.loc[by_segment.n.idxmax(), 'miles']:.3f} miles")
print(f"sliding 1.0 mile : max {window.max()} structures, "
      f"starting at postmile {starts[peak]:.1f}")
print(by_segment.to_string(index=False))

fixed miles      :  70 bins, max 10 structures
route segments   :   9 bins, max 26 structures in 17.979 miles
sliding 1.0 mile : max 10 structures, starting at postmile 15.9
 from_pm  to_pm  miles  n  deck_ft
   0.000  8.108  8.108 10   1001.5
   8.198 19.538 11.340 22   1071.5
  19.548 25.817  6.269  0      0.0
  25.825 32.976  7.151 17    728.2
  32.976 34.567  1.591  0      0.0
  34.616 49.063 14.447 15    951.6
  49.067 49.471  0.404  3    150.5
  49.477 67.456 17.979 26   2222.2
  67.456 69.322  1.866  0      0.0

Three segments come back empty, and all three are the R-prefixed ones. The bridge inventory's PM column carries no prefix, so its postmiles were matched against ranges they do not belong to — a reminder that binning by segment inherits whatever key problem the event table has. The fixed-mile binning hides that by not using the key at all, which is not the same as being right.

Example 2 — weighting by what the events carry

per_mile = (on101.assign(mile=np.floor(on101.PM).astype(int))
            .groupby("mile").agg(n=("PM", "size"), deck_ft=("LENG", "sum")))
per_mile = per_mile.reindex(range(70), fill_value=0)

print("by count:")
print(per_mile.nlargest(4, "n").to_string())
print("\nby deck length:")
print(per_mile.nlargest(4, "deck_ft").to_string())
overlap = set(per_mile.nlargest(5, "n").index) & set(per_mile.nlargest(5, "deck_ft").index)
print(f"\nmiles in both top fives: {sorted(overlap)}")

by count:
       n  deck_ft
mile             
16    10    570.0
28     5    179.5
29     5    237.2
45     5    158.2

by deck length:
       n  deck_ft
mile             
55     3    819.6
0      3    705.8
16    10    570.0
35     2    405.4

miles in both top fives: [16]

Example 3 — a rate per vehicle-mile

aadt = pd.DataFrame([f["attributes"] for f in
                     query("Traffic_AADT", "CNTY='SLO' AND RTE='101'")])
aadt["PM"] = pd.to_numeric(aadt.PM, errors="coerce")
aadt["AHEAD_AADT"] = pd.to_numeric(aadt.AHEAD_AADT, errors="coerce")
st = (aadt.drop_duplicates(subset=["PM_PFX", "PM", "DESCRIPTION"])
      .sort_values("PM").dropna(subset=["PM"]).reset_index(drop=True))

def volume_at(pm):
    prior = st[st.PM <= pm]
    return float(prior.AHEAD_AADT.iloc[-1]) if len(prior) else np.nan

rows = []
for mile in range(70):
    vol = volume_at(mile + 0.5)
    n = int(per_mile.n.get(mile, 0))
    vmt = vol * 1.0 * 365 / 1e6 if np.isfinite(vol) else np.nan     # million vehicle-miles/yr
    rows.append({"mile": mile, "structures": n, "aadt": vol,
                 "per_100m_vmt": 100 * n / vmt if vmt else np.nan})
rate = pd.DataFrame(rows).dropna(subset=["aadt"])

print(rate.nlargest(4, "structures")[["mile", "structures", "aadt", "per_100m_vmt"]]
      .round(2).to_string(index=False))
print()
print(rate.nlargest(4, "per_100m_vmt")[["mile", "structures", "aadt", "per_100m_vmt"]]
      .round(2).to_string(index=False))

 mile  structures    aadt  per_100m_vmt
   16          10 71000.0         38.59
   28           5 65000.0         21.07
   29           5 49000.0         27.96
   45           5 56000.0         24.46

 mile  structures    aadt  per_100m_vmt
   65           4 20000.0         54.79
   16          10 71000.0         38.59
   29           5 49000.0         27.96
   67           2 20000.0         27.40

The two rankings disagree, which is the point of normalising: the mile with the most structures is a busy urban stretch, and the mile with the most structures per vehicle-mile is a quiet one where the same count is spread over under a third of the traffic.

Why sliding windows find more

Why the sliding window finds more

A fixed grid puts a boundary every mile whether or not anything is happening there. A cluster that straddles one of those boundaries is split between two bins and neither looks remarkable. Sliding the window in small steps means every possible position is tested, so the maximum is the real maximum — at the cost of overlapping windows that cannot be summed.

Why counts need a denominator

Counts per mile answer "where are the most of these", which is usually a proxy for "where is the most road" or "where is the most traffic". Dividing by exposure — length, vehicle-miles, population — turns it into a rate that can be compared between places. Which denominator depends on the question, and the wrong one produces a map of traffic volume with a different legend.

Why the bin edges belong in the output

Two analyses of the same events with different bins produce different peaks, different maxima and different conclusions, and nothing in either result says which bin was used. Half-mile bins, one-mile bins and the route's own segmentation all appear above; only the third is reproducible without being told the parameters.

Three binnings of the same 93 structures, and three different busiest stretches.

Boundary events, zero exposure and small numbers

Events on the boundary between bins: pick half-open intervals and say which end is open.
Linear events need apportioning by overlap, not by their midpoint, when they cross a bin.
Zero-exposure bins (no road, or no volume data) should be dropped rather than counted as zero.
Small numbers: a rate from two events in a mile has an enormous confidence interval; say so.
Both carriageways: decide whether a divided highway's structures are one route or two before counting.
Seasonal volumes: AADT is an annual average, so a rate built on it is an annual rate.

Internal links

Dynamic segmentation explained — where segment boundaries come from
How to overlay linear event tables in Python — combining the tables first
How to find the measure of a point on a route — getting events onto the route
How to aggregate spatial data by region in GeoPandas — the polygon version
Event lengths do not add up to the route length — checking the denominator
Hazard, exposure and vulnerability explained — the same normalisation argument

FAQ

How do I count events per mile on a route?

Bin the measures with np.histogram or a groupby on the floor of the measure, and state the bin width.

Fixed bins or the route's own segments?

Fixed bins for reporting, route segments for joining to other attributes; they give different answers.

What is a sliding window good for?

Finding concentrations that a fixed grid would split. Its windows overlap, so the counts cannot be summed.

Should I divide by length or by traffic?

By length to compare stretches of road, by vehicle-miles to compare exposure. Say which.

What about bins with no road?

Drop them. A gap in the numbering is not a stretch with zero events.

How do I weight by event size?

Sum an attribute — deck area, cost, lane-miles — instead of counting rows.