How to extract the part of a route between two measures

Problem statement

Line events are stored as a route and two measures: a speed limit from postmile 12.4 to 18.9, a resurfacing contract from 30.1 to 37.9, a traffic volume between two count stations. Drawing them means cutting the route between those measures — and the cut has to respect the calibration, cross segment boundaries where the event does, and report honestly when part of the range has no route under it.

This guide turns 47 traffic volume events on US 101 into geometry, and checks that the pieces add up.

Quick answer

Convert both measures to distances, then take the substring:

import numpy as np
from shapely.ops import substring
from lrs import routes

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

def along(line):
    xy = np.asarray(line.coords)
    return np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))])

def cut(seg, m0, m1):
    """The piece of one segment between two measures, in its own direction."""
    m, d = np.asarray(seg.m), along(seg.geometry)
    order = np.argsort(m)
    d0, d1 = (float(np.interp(v, m[order], d[order])) for v in (m0, m1))
    return substring(seg.geometry, min(d0, d1), max(d0, d1))

seg = shn.iloc[0]
piece = cut(seg, 2.0, 6.0)
print(f"segment {seg.bPM:.3f}–{seg.ePM:.3f}, {seg.geometry.length / 1609.344:.3f} mi")
print(f"postmiles 2.0–6.0 -> {piece.length / 1609.344:.4f} miles of geometry, "
      f"{len(piece.coords)} vertices")
print(f"  measure span requested: 4.000 mi; ratio {piece.length / 1609.344 / 4:.4f}")

segment 0.000–8.108, 8.118 mi
postmiles 2.0–6.0 -> 3.9959 miles of geometry, 41 vertices
  measure span requested: 4.000 mi; ratio 0.9990

The extracted length is not exactly four miles, and it should not be: four miles of measure is whatever length of road the calibration says lies between those two numbers.

Clip to each segment, convert to distances, cut separately

1. Clip the request to the segment

An event may start before a segment and end after it. Intersect the two ranges first and cut only what overlaps.

2. Convert measures to distances, not fractions

substring takes distances along the line. Getting them from np.interp against the stored measures is what keeps the cut calibrated.

3. Cut each overlapping segment separately

An event spanning three segments produces three pieces. Keep them as separate rows or join them into a MultiLineString; do not pretend it is one line.

4. Account for the measure range that had no route

Sum the covered measure and compare it with the requested span. The difference is gaps, prefix mismatches or measures off the end.

5. Preserve direction if it matters

substring returns the piece in the line's own direction. For drawing it makes no difference; for anything that reads the first or last vertex it does.

6. Check the total

Cutting every event and summing the lengths should reproduce the route's length. If it does not, the events overlap or leave gaps.

Cutting events that cross segments

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 — cutting an event that crosses segments

def cut_event(frame, m0, m1, prefix=""):
    """Every geometry piece for a measure range, with the measure each covers."""
    lo, hi = sorted((m0, m1))
    pieces = []
    for seg in frame.itertuples():
        if (seg.PMPrefix or "") != prefix:
            continue
        m = np.asarray(seg.m)
        s_lo, s_hi = min(m[0], m[-1]), max(m[0], m[-1])
        a, b = max(lo, s_lo), min(hi, s_hi)
        if b - a > 1e-9:
            pieces.append({"from_pm": a, "to_pm": b, "measure_mi": b - a,
                           "geometry": cut(seg, a, b)})
    return pieces

for span in [(2.0, 6.0), (7.0, 10.0), (18.0, 27.0)]:
    pieces = cut_event(shn, *span)
    measure = sum(p["measure_mi"] for p in pieces)
    length = sum(p["geometry"].length for p in pieces) / 1609.344
    print(f"postmiles {span[0]}–{span[1]}: {len(pieces)} piece(s), "
          f"{measure:.3f} mi of measure covered of {span[1] - span[0]:.3f} requested, "
          f"{length:.3f} mi of geometry")

postmiles 2.0–6.0: 1 piece(s), 4.000 mi of measure covered of 4.000 requested, 3.996 mi of geometry
postmiles 7.0–10.0: 2 piece(s), 2.910 mi of measure covered of 3.000 requested, 2.908 mi of geometry
postmiles 18.0–27.0: 2 piece(s), 2.713 mi of measure covered of 9.000 requested, 2.713 mi of geometry

The 18.0–27.0 request covers 2.713 miles of the nine asked for: postmiles 19.548 to 25.817 are under the R prefix, which this call did not ask for, so only the unprefixed ends of the range come back. The covered-measure figure is what surfaces that; the geometry on its own looks like a perfectly reasonable pair of short lines.

Example 2 — a whole event table as geometry

import geopandas as gpd, pandas as pd
from lrs import query

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")
aadt["prefix"] = aadt.PM_PFX.fillna("")
stations = (aadt.drop_duplicates(subset=["prefix", "PM", "DESCRIPTION"])
            .sort_values("PM").reset_index(drop=True))

rows = []
for a, b in zip(stations.itertuples(), stations.shift(-1).itertuples()):
    if pd.isna(b.PM) or b.PM <= a.PM or pd.isna(a.AHEAD_AADT):
        continue
    for piece in cut_event(shn, a.PM, b.PM, a.prefix):
        rows.append({**piece, "aadt": a.AHEAD_AADT, "from_station": a.DESCRIPTION})

events = gpd.GeoDataFrame(rows, crs=shn.crs)
print(f"{len(stations) - 1} volume events -> {len(events)} geometry pieces")
print(f"measure covered {events.measure_mi.sum():.3f} mi of "
      f"{stations.PM.max() - stations.PM.min():.3f}")
print(f"geometry length  {events.geometry.length.sum() / 1609.344:.3f} mi, "
      f"route length {shn.geometry.length.sum() / 1609.344:.3f} mi")
print(events.nlargest(3, "aadt")[["from_pm", "to_pm", "aadt", "from_station"]]
      .to_string(index=False))

47 volume events -> 51 geometry pieces
measure covered 66.958 mi of 69.322
geometry length  66.966 mi, route length 69.171 mi
 from_pm  to_pm    aadt                              from_station
  15.579 16.398 84000.0                   PISMO BEACH, PISMO OAKS
  14.613 15.579 75000.0                PISMO BEACH, OAK PARK ROAD
   0.000  0.813 74000.0 SANTA BARBARA/SAN LUIS OBISPO COUNTY LINE

Example 3 — do the pieces add up?

overlap = 0.0
ordered = events.sort_values("from_pm")
prev_to = None
for r in ordered.itertuples():
    if prev_to is not None and r.from_pm < prev_to - 1e-9:
        overlap += prev_to - r.from_pm
    prev_to = max(prev_to or r.to_pm, r.to_pm)

covered = events.measure_mi.sum()
span = ordered.to_pm.max() - ordered.from_pm.min()
print(f"events span postmiles {ordered.from_pm.min():.3f}–{ordered.to_pm.max():.3f} "
      f"({span:.3f} mi)")
print(f"measure covered by pieces: {covered:.3f} mi")
print(f"overlap between pieces:    {overlap:.3f} mi")
print(f"uncovered inside the span: {span - covered + overlap:.3f} mi")

events span postmiles 0.000–69.322 (69.322 mi)
measure covered by pieces: 66.958 mi
overlap between pieces:    0.000 mi
uncovered inside the span: 2.364 mi

Two and a third miles uncovered, and no overlap at all. The uncovered measure is spread over seven events, and the largest single piece is 1.640 miles: the event from postmile 30.360 to 37.863 is unprefixed and the road between 32.976 and 34.567 is numbered R, so the middle of that event has no unprefixed route under it. The rest are the equation gaps — 0.090 of a mile at 8.108, a few thousandths at each of the others.

Every line-event conversion should print these three numbers. They are the difference between a map that looks complete and one that is, and here they say that 3% of a route's volume events cannot be drawn against the route they name.

Substring rather than split

Why substring and not split

split cuts at a geometry and returns the pieces on either side, which needs a point exactly on the line and returns two things when you wanted one. substring(line, start, end) takes distances and returns the piece between them, interpolating new end vertices where needed — which is the operation this is.

Why measure length and geometry length differ

The measure span is a number from the register; the geometry length is a property of the digitised line. They agree to within a percent on most segments here and not exactly on any, because the two were produced by different processes decades apart. Report whichever the question asks for, and say which.

Why the uncovered measure is the thing to watch

A cut that returns geometry always produces a map. The map has no way of showing measure that had no route under it — the missing tenth of a mile simply is not drawn, and the line either side looks continuous. The only way to see it is to compare covered measure with requested measure, which is why every function here returns both.

Reversed measures and very short events

• Reversed measures: sort the measure and distance arrays together; substring needs ascending distances.
• Zero-length results happen when both measures fall between the same two vertices; keep or drop them deliberately.
• normalized=True on substring takes fractions, which is exactly the proportional assumption to avoid.
• MultiLineString results from joining pieces are fine for drawing and awkward for length-weighted joins.
• Very short events — a bridge deck, a gore area — may be shorter than the vertex spacing; the result is a two-point line.
• PostGIS ST_LocateBetween cuts on the M ordinate directly and returns a collection.

Internal links

• How to find the point at a given measure on a route — the point form of this
• Dynamic segmentation explained — where the measure ranges come from
• How to overlay linear event tables in Python — combining several of them
• Event lengths do not add up to the route length — the fix guide
• Route calibration explained — why four miles of measure is not four miles of line
• How to clean a line network: dangles, overshoots and disconnects — the geometry side

FAQ

How do I cut a line between two measures?

Interpolate both measures to distances along the line, then call shapely.ops.substring with them.

What if the event spans several segments?

Cut each overlapping segment separately and keep the pieces, along with the measure each one covers.

Why is the piece not the length I asked for?

Because the measure span and the geometry length are different quantities on a calibrated route.

What happens to measures with no route under them?

Nothing is drawn. Compare covered measure with requested measure to find them.

Can I use substring with fractions?

You can, and it assumes measures are proportional to length, which is the error this whole approach avoids.

How do I check the result?

Sum the covered measure, the overlap and the requested span; they should reconcile exactly.