How to overlay linear event tables in Python
On this page
Problem statement
Three tables describe the same stretch of highway: traffic volumes between count stations, a functional classification published against a different measure system, and the route's own prefix segments. Answering "what is the volume on the principal-arterial parts of this route?" means overlaying them.
The overlay itself is twenty lines. Everything around it — translating one table's measures, keeping the route key, deciding what a piece covered by only two of the three tables means — is where the work is.
Quick answer
One function: collect breakpoints, cut, label each piece from every table by midpoint.
import numpy as np, pandas as pd
def overlay(tables, start, end, tol=1e-6):
"""tables: {name: [(from_m, to_m, value), …]}. Returns homogeneous pieces."""
cuts = {round(start, 6), round(end, 6)}
for rows in tables.values():
for a, b, _ in rows:
if end > a and b > start:
cuts.add(round(max(a, start), 6))
cuts.add(round(min(b, end), 6))
edges = sorted(cuts)
pieces = []
for lo, hi in zip(edges[:-1], edges[1:]):
if hi - lo <= tol:
continue
mid = (lo + hi) / 2
row = {"from_m": lo, "to_m": hi, "length": hi - lo}
for name, rows in tables.items():
row[name] = next((v for a, b, v in rows if a <= mid < b), None)
pieces.append(row)
return pd.DataFrame(pieces)
demo = overlay({"speed": [(0, 5, 55), (5, 12, 65)],
"surface": [(0, 3, "asphalt"), (3, 12, "concrete")]}, 0, 12)
print(demo.to_string(index=False)) from_m to_m length speed surface
0 3 3 55 asphalt
3 5 2 55 concrete
5 12 7 65 concrete
Two tables with one breakpoint each produce three pieces. That is the whole idea; the rest of this guide is about making it survive real data.
One measure system, clipped spans, rounded breakpoints, midpoints
1. Get every table onto one measure system
Tables published against different measures have to be translated before the breakpoints mean anything. A crosswalk from the route table is the safe way.
2. Clip to the route span
An event table may cover more than the route you are analysing. Clipping at the start makes the accounting at the end interpretable.
3. Union the breakpoints, with rounding
Measures are published to a fixed precision. Round the cuts to that precision or floating point will produce pieces a micron long.
4. Label by midpoint
a <= mid < b on the midpoint of each piece. No boundary case, no half-open-interval argument.
5. Leave missing values missing
A piece no table covers gets None, not a filled-forward value. The count of those is the quality measure.
6. Aggregate with length weights
Once the pieces exist, every summary is a groupby with length as the weight.
Three real tables on one route
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 real tables on one route
from lrs import query, routes
shn = routes("County='SLO' AND Route=101 AND AlignCode='Right'").sort_values("bPM")
cross = shn[["PMPrefix", "bPM", "ePM", "bOdometer", "eOdometer"]].copy()
cross["PMPrefix"] = cross.PMPrefix.fillna("")
def odometer_to_pm(odo):
for r in cross.itertuples():
if r.bOdometer <= odo <= r.eOdometer and r.eOdometer > r.bOdometer:
frac = (odo - r.bOdometer) / (r.eOdometer - r.bOdometer)
return r.bPM + frac * (r.ePM - r.bPM)
return None
# volumes, keyed on postmile already
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))
volume = [(a.PM, b.PM, a.AHEAD_AADT)
for a, b in zip(st.itertuples(), st.shift(-1).itertuples())
if not pd.isna(b.PM) and b.PM > a.PM and not pd.isna(a.AHEAD_AADT)]
# functional class, keyed on the odometer: translate it
fclass_raw = []
for f in query("CRS_Functional_Classification", "RouteID='SHS_101._P'", want_m=True):
for path in f["geometry"]["paths"]:
fclass_raw.append((path[0][2], path[-1][2], f["attributes"]["F_System"]))
fclass = []
for a, b, v in fclass_raw:
pa, pb = odometer_to_pm(min(a, b)), odometer_to_pm(max(a, b))
if pa is not None and pb is not None and pb > pa:
fclass.append((pa, pb, v))
prefix = [(s.bPM, s.ePM, (s.PMPrefix or "(none)")) for s in shn.itertuples()]
print(f"volume events {len(volume)}, functional-class pieces "
f"{len(fclass)} of {len(fclass_raw)} translated, route segments {len(prefix)}")
pieces = overlay({"aadt": volume, "fclass": fclass, "prefix": prefix},
shn.bPM.min(), shn.ePM.max())
print(f"-> {len(pieces)} homogeneous pieces covering "
f"{pieces.length.sum():.3f} miles")
print(pieces.head(6).to_string(index=False))volume events 47, functional-class pieces 164 of 2285 translated, route segments 9
-> 226 homogeneous pieces covering 69.322 miles
from_m to_m length aadt fclass prefix
0.000000 0.537963 0.537963 74000.0 NaN (none)
0.537963 0.767961 0.229998 74000.0 2.0 (none)
0.767961 0.813000 0.045039 74000.0 2.0 (none)
0.813000 1.014958 0.201958 65000.0 2.0 (none)
1.014958 4.637377 3.622419 65000.0 2.0 (none)
4.637377 4.807410 0.170033 65000.0 2.0 (none)
164 of the 2,285 functional-class pieces translate: the rest belong to the other thirteen counties on this route, whose odometer values fall outside the crosswalk built from San Luis Obispo's segments. Translating first and printing the survivors is what makes that a number rather than a surprise.
Example 2 — what the overlay tells you that the tables do not
have = pieces.notna()
for col in ("aadt", "fclass", "prefix"):
covered = pieces.loc[have[col], "length"].sum()
print(f"{col:8} covers {covered:7.3f} mi of {pieces.length.sum():.3f} "
f"({covered / pieces.length.sum():5.1%}) in {int(have[col].sum()):4} pieces")
complete = pieces[have[["aadt", "fclass", "prefix"]].all(axis=1)]
print(f"\npieces with all three: {len(complete)}, {complete.length.sum():.3f} mi")
partial = pieces[~have[["aadt", "fclass", "prefix"]].all(axis=1)]
print(f"pieces missing something: {len(partial)}, {partial.length.sum():.3f} mi")
print(partial.head(5).to_string(index=False))aadt covers 69.322 mi of 69.322 (100.0%) in 226 pieces
fclass covers 68.753 mi of 69.322 (99.2%) in 224 pieces
prefix covers 69.155 mi of 69.322 (99.8%) in 220 pieces
pieces with all three: 218, 68.586 mi
pieces missing something: 8, 0.736 mi
from_m to_m length aadt fclass prefix
0.000 0.537963 0.537963 74000.0 NaN (none)
8.108 8.198000 0.090000 64000.0 2.0 NaN
19.538 19.548000 0.010000 71000.0 2.0 NaN
25.817 25.825000 0.008000 69000.0 2.0 NaN
34.567 34.616000 0.049000 48000.0 2.0 NaN
Example 3 — the question the overlay was for
by_class = (complete.assign(vmt=complete.aadt * complete.length)
.groupby("fclass")
.agg(miles=("length", "sum"), vmt=("vmt", "sum"),
pieces=("length", "size")))
by_class["mean_aadt"] = by_class.vmt / by_class.miles
print(by_class.round(1).to_string())
by_prefix = (complete.assign(vmt=complete.aadt * complete.length)
.groupby("prefix").agg(miles=("length", "sum"), vmt=("vmt", "sum")))
by_prefix["mean_aadt"] = (by_prefix.vmt / by_prefix.miles).round(0)
print("\n" + by_prefix.round(1).to_string()) miles vmt pieces mean_aadt
fclass
2.0 68.6 3711037.5 218 54107.9
miles vmt mean_aadt
prefix
(none) 58.9 3154631.7 53567.0
R 9.7 556405.8 57392.0
Length-weighted means, computed from pieces that are homogeneous in all three attributes. No piece contributes to two classes, and no class borrows a volume from a stretch that has a different one.
The functional-class table turns out to have a single value along this whole route in this county — every piece is class 2, principal arterial — so that summary has one row. That is a result, not a wasted overlay: the alternative was assuming it, and the prefix split underneath shows the realigned nine miles carrying slightly more traffic than the rest.
Translate first, and leave gaps empty
Why the translation belongs outside the overlay
A table in another measure system is not an overlay problem, it is a conflation problem, and mixing the two hides where events were lost. Translating first, counting what survived, then overlaying means the overlay's coverage numbers describe the overlay rather than the crosswalk.
Why None is better than a fill
Filling a missing attribute forward is a modelling decision dressed as a data-cleaning step. On a route where the classification stops at a county line, a forward fill extends it silently into a county that classified it differently. Leaving the gap means the summary explicitly excludes it, and the excluded mileage is printed.
Why this scales
The overlay is O(pieces × tables) with a linear scan per lookup, which is fine for a route and slow for a state. For a large run, sort each table once and walk the breakpoints with an index per table, or push the whole thing into a database with range joins — the logic is the same and the midpoint rule still does the labelling.
Duplicates, overlaps and dated rows
- Duplicate rows in an event table produce duplicate breakpoints, which the rounding-plus-set handles.
- Overlapping rows within one table make the midpoint lookup return the first match; detect them first.
- Direction-specific tables need the alignment in the key, or they will overlay onto both carriageways.
- Dates: filter temporal events to one point in time before overlaying, or a work zone from 2019 will cut the route.
- Very many tables make the piece count multiply; overlay only the tables the question needs.
- SQL does this with range joins;
BETWEENon the midpoint is the same rule.
Internal links
- Dynamic segmentation explained — what this operation is
- LRS conflation explained — translating a table into the right measures
- How to extract the part of a route between two measures — turning the pieces into geometry
- How to aggregate events to route segments — the summary step
- Event lengths do not add up to the route length — reading the coverage numbers
- How to aggregate spatial data by region in GeoPandas — the polygon equivalent
FAQ
How do I combine two linear event tables?
Union their breakpoints, cut the route at each one, and label every piece by looking up its midpoint in each table.
What if the tables use different measures?
Translate one into the other's system first, and count how many events survived the translation.
Why use the midpoint?
It is unambiguously inside exactly one interval of each table, so there is no boundary case.
What should a piece with no value contain?
None. Filling it in invents data, and the uncovered mileage is the number worth reporting.
How do I summarise the result?
Group by the attributes you care about and weight by piece length.
Does this work for point events?
Point events do not cut the route; aggregate them onto the pieces afterwards.