Route measures are not monotonic

Problem statement

np.interp returns nonsense, events land at the wrong end of the route, and a plot of measure against distance along the line looks like a saw. The measures are not increasing along the geometry, and every operation that assumes they are — interpolation, substring, binning — is quietly wrong rather than loudly broken.

There are four causes, and they need different fixes: the pieces were assembled in the wrong order, the geometry is digitised against the numbering, the register genuinely repeats measures, or consecutive vertices are duplicates.

Quick answer

Classify every path before using it, and look at what the classes are:

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

shn = pd.concat([routes(f"District={d}") for d in (4, 5, 7)], ignore_index=True)

def classify(m):
    d = np.diff(np.asarray(m))
    if np.all(d > 0):
        return "strictly increasing"
    if np.all(d >= 0):
        return "increasing with ties"
    if np.all(d <= 0):
        return "decreasing"
    return "not monotonic"

shn["shape"] = [classify(m) for m in shn.m]
print(f"{len(shn)} paths")
print(shn["shape"].value_counts().to_string())

1633 paths
shape
strictly increasing     1271
increasing with ties     362

A published LRS is usually clean, and checking costs one line. The measures that go wrong are almost always the ones you built yourself.

Plot measure against distance, reverse rather than sort

1. Plot measure against distance along the line

The shape names the cause: a saw is the wrong assembly order, a downward slope is reversed geometry, a step is a repeated measure, a flat spot is duplicate vertices.

2. Reverse, do not sort, a decreasing path

A decreasing path is correct data with the geometry running the other way. Reverse both arrays together; sorting the measures alone destroys the correspondence with the vertices.

3. Re-order assembled pieces by connectivity

If you merged the route yourself, the fix is upstream: order the pieces head-to-tail before concatenating, and reverse any piece whose start is nearer the previous piece's start than its end.

4. Drop duplicate vertices

Consecutive identical coordinates produce a zero-length step and a tied measure. Removing them is safe; it changes no geometry.

5. Split at a genuine repeat

If the register really does use a measure twice — an overlap in the numbering — the path is two routes, not one. Split it and give the pieces their prefixes.

6. Re-check after every repair

The repairs interact: reversing a merged route can turn a saw into a different saw. Classify again and only then interpolate.

A sawtooth built by accident, and its repair

The data used on this page

The street network is the OpenStreetMap driving network for inner Portland, Oregon, built with OSMnx and saved as GraphML. overpass-api.de refuses connections from here, so the query went to the maps.mail.ru mirror of the same API, and the setup reads that answer as saved on 9 October 2026: OpenStreetMap is edited continuously, and a live query moves the edge counts a little.

import io
import os
import requests
import zipfile

import osmnx as ox

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-measures-not-monotonic-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # overpass-api.de mirror
ox.settings.overpass_rate_limit = False
PDX = (-122.75, 45.46, -122.56, 45.57)                 # west, south, east, north: inner Portland
if not os.path.exists("pdx_drive.graphml"):
    ox.save_graphml(ox.graph_from_bbox(PDX, network_type="drive"), "pdx_drive.graphml")
G = ox.load_graphml("pdx_drive.graphml")
print(f"drive network: {len(G):,} nodes, {G.number_of_edges():,} edges")

drive network: 13,757 nodes, 39,806 edges

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 — a real sawtooth, built by accident

import osmnx as ox
from shapely.ops import linemerge

edges = ox.graph_to_gdfs(ox.load_graphml("pdx_drive.graphml"), nodes=False).to_crs(3857)
street = edges[edges.name.astype(str) == "Southeast Division Street"]

seen, unique = set(), []
for g in street.geometry:
    key = tuple(sorted([tuple(np.round(g.coords[0], 2)), tuple(np.round(g.coords[-1], 2))]))
    if key not in seen:
        seen.add(key)
        unique.append(g)

# the wrong way: concatenate the pieces in the order they happen to arrive
naive_xy = np.vstack([np.asarray(g.coords) for g in unique])
naive_m = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(naive_xy, axis=0).T))]) / 1000

steps = np.diff(naive_m)
jump = np.hypot(*np.diff(naive_xy, axis=0).T)
print(f"naive concatenation: {len(naive_xy)} vertices, 'route' {naive_m[-1]:.3f} km long")
print(f"  steps over 200 m between consecutive vertices: {(jump > 200).sum()}")
print(f"  largest jump: {jump.max():,.0f} m — the gap between two unconnected pieces")

merged = linemerge(unique)
parts = sorted(merged.geoms, key=lambda g: -g.length) if merged.geom_type == "MultiLineString" else [merged]
route = parts[0]
good_xy = np.asarray(route.coords)
good_m = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(good_xy, axis=0).T))]) / 1000
print(f"after linemerge: {len(good_xy)} vertices, {good_m[-1]:.3f} km, "
      f"largest step {np.hypot(*np.diff(good_xy, axis=0).T).max():.0f} m")

naive concatenation: 574 vertices, 'route' 281.201 km long
  steps over 200 m between consecutive vertices: 71
  largest jump: 10,083 m — the gap between two unconnected pieces
after linemerge: 294 vertices, 8.514 km, largest step 115 m

The naive measures are monotonic — cumulative length always is — and completely wrong: the "route" is 281 km long, with 71 jumps of more than 200 m between unconnected pieces, one of them ten kilometres. Monotonic is necessary, not sufficient, and the test that catches this is the step between consecutive vertices rather than the direction of the measures.

Example 2 — reversing rather than sorting

fake_m = good_m[::-1].copy()                    # as if the geometry ran the other way
print(f"decreasing measures: {fake_m[0]:.3f} -> {fake_m[-1]:.3f}")

probe = 2.0
wrong = float(np.interp(probe, fake_m, np.arange(len(fake_m))))
xy_rev, m_rev = good_xy[::-1], fake_m[::-1]
right = float(np.interp(probe, m_rev, np.arange(len(m_rev))))
print(f"np.interp on decreasing x: vertex {wrong:.0f} (np.interp assumed increasing)")
print(f"after reversing both arrays: vertex {right:.0f} of {len(fake_m)}")

sorted_only = np.sort(fake_m)
print(f"sorting the measures alone: first three {sorted_only[:3].round(4).tolist()}, "
      f"but they no longer match the vertices")

decreasing measures: 8.514 -> 0.000
np.interp on decreasing x: vertex 293 (np.interp assumed increasing)
after reversing both arrays: vertex 69 of 294
sorting the measures alone: first three [0.0, 0.0412, 0.0647], but they no longer match the vertices

np.interp does not raise on decreasing x. It returns a number, and the number is wrong in a way that looks like an ordinary measure.

Example 3 — ties and duplicate vertices

ties = []
for r in shn.itertuples():
    m = np.asarray(r.m)
    d = np.diff(m)
    if np.any(d == 0):
        xy = np.asarray(r.geometry.coords)
        step = np.hypot(*np.diff(xy, axis=0).T)
        ties.append({"county": r.County, "route": r.Route,
                     "tied_steps": int((d == 0).sum()),
                     "flat_metres": float(step[d == 0].sum())})
ties = pd.DataFrame(ties)
print(f"{len(ties)} of {len(shn)} paths have at least one tied measure")
print(f"total geometry under a tied measure: {ties.flat_metres.sum():.1f} m; "
      f"worst single path {ties.flat_metres.max():.1f} m")
print(ties.nlargest(3, "flat_metres").to_string(index=False))

362 of 1633 paths have at least one tied measure
total geometry under a tied measure: 29.6 m; worst single path 0.7 m
county  route  tied_steps  flat_metres
   SON    128          10     0.735065
   SON    128          10     0.735065
   SON     12          12     0.514196

Thirty metres across three districts, and the worst single path has 0.7 m of it: these are duplicate vertices, not a numbering problem, and dropping them is safe. A tie covering hundreds of metres would be something else — a stretch of road where the register stops advancing, which cannot be interpolated through at all.

Why np.interp stays silent

Why np.interp will not tell you

np.interp documents that x must be increasing and does not check. With decreasing x it returns the value at the first bracketing pair it finds, which for a fully reversed array is the endpoint. No exception, no warning, and an answer in the right range — the worst possible failure mode for a function used a million times in a batch job.

Why sorting is the wrong repair

The measures and the vertices are two halves of one structure. Sorting the measures alone silently re-pairs measure i with vertex j, which produces a monotonic array describing a route that does not exist. Every repair has to move both arrays, or neither.

Why assembled routes are the usual culprit

A published LRS is maintained by people whose job is the register, and it is nearly always monotonic — in 1,633 paths across three Caltrans districts, every one increases. A route assembled from OSM ways, a shapefile of fragments, or a spatial join is assembled by ordering rules that break at every junction, and the measures inherit the breakage.

Near-ties, loops and multipart routes

• Floating point makes exact ties rare and near-ties common; compare with a tolerance the register's precision justifies.
• Loops and switchbacks are monotonic in measure and cross themselves in space; that is fine.
• MultiLineString routes have no single vertex order; explode them and measure each part.
• Measures in a joined table can be non-monotonic even when the geometry is fine — check both.
• np.interp with period= is for angles, not for this.
• Check after every geometry edit: simplify, snap and densify can all reorder or drop vertices.

Internal links

• M values explained: geometry that carries a measure — what the arrays are
• How to build measures onto a line that has none — where sawtooth measures come from
• Route calibration explained — what monotonic measures are supposed to mean
• Event lengths do not add up to the route length — the other arithmetic surprise
• How to clean a line network: dangles, overshoots and disconnects — fixing the geometry first
• How to read a linear referencing system in Python — the check as part of loading

FAQ

Why are my route measures not increasing?

The pieces were assembled out of order, the geometry runs against the numbering, the register repeats a measure, or there are duplicate vertices.

Can I just sort the measures?

No. Sorting breaks the pairing between measures and vertices. Reverse both arrays, or reorder both.

Does np.interp warn about decreasing x?

No. It returns a plausible wrong number.

Are tied measures a problem?

Only if they cover real distance. Duplicate vertices produce harmless ties a few centimetres long.

How do I know a published LRS is fine?

Classify every path when you load it; it costs one line and it is the cheapest quality check available.

What if the register really repeats a measure?

Then the path is two routes. Split it and carry the prefix or suffix that distinguishes them.