⚙️ How-To Linear referencing Networks & routing Linear Referencing · step 9 Notebook By · Updated

How to find the point at a given measure on a route

On this page

Problem statement

An event table gives a route and a measure and no coordinates — which is the normal state of highway asset data. Putting it on a map means the inverse of locating a point: find the segment whose measure range contains the value, find the distance along that segment where the measure reaches it, and interpolate the position.

Each of those three steps has a failure mode. The wrong segment is chosen because the prefix was ignored; the measure falls in a gap and there is no segment at all; and the distance is computed from proportion instead of from the calibration.

Quick answer

Select on the whole key, then interpolate distance from measure:

import numpy as np
from lrs import routes

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

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

def point_at(frame, pm, prefix=""):
    """The position at a measure, or None if no segment carries it."""
    for seg in frame.itertuples():
        m = np.asarray(seg.m)
        lo, hi = min(m[0], m[-1]), max(m[0], m[-1])
        if (seg.PMPrefix or "") == prefix and lo - 1e-9 <= pm <= hi + 1e-9:
            d = along_distances(seg.geometry)
            order = np.argsort(m)                       # np.interp needs increasing x
            return seg.geometry.interpolate(float(np.interp(pm, m[order], d[order])))
    return None

for pm, prefix in [(49.64, ""), (8.15, ""), (20.0, ""), (20.0, "R")]:
    p = point_at(shn, pm, prefix)
    label = f"{prefix}{pm}"
    print(f"postmile {label:>7}: " +
          ("no segment carries it" if p is None else f"({p.x:,.0f}, {p.y:,.0f})"))
postmile   49.64: (-64,127, -276,072)
postmile    8.15: no segment carries it
postmile    20.0: no segment carries it
postmile   R20.0: (-62,722, -316,519)

Two of the four have no answer, and returning None is the right one. A function that falls back to the nearest segment would put an event a mile from where its record says it is.

Map of US-101 through San Luis Obispo County with postmile ticks every five miles; most of the road has plain postmiles, three short sections carry the R prefix, and red dots mark postmiles 49.64 and R20.0. Map of US-101 through San Luis Obispo County with postmile ticks every five miles; most of the road has plain postmiles, three short sections carry the R prefix, and red dots mark postmiles 49.64 and R20.0.
The route the code queries, with a tick every five postmiles. The two red dots are the answers. PM 8.15 has none because it falls in the 0.09-mile equation gap between segments ending at 8.108 and starting at 8.198; plain PM 20.0 has none because that stretch was realigned and is numbered R19.548 to R25.817. Only the R-prefixed query finds it.
From a route key and a measure to a coordinate. measure to coordinate1select the segment on the full keythe whole key, not the number2sort measures and distances togethernp.interp needs increasing x3interpolate the distance alongnp.interp(pm, measures, distances)4interpolate the positionline.interpolate(distance)5return None where nothing carries ita gap has no positionFour of these five lines are lookups; the first is where it goes wrong silently spatialworkflow.io
Three lookups, and the first one is the one that goes wrong silently.

Full key, sorted measures, interpolate, return nothing in gaps

1. Select the segment on the full key

County, route, prefix, suffix and alignment. A measure that exists under two prefixes needs the prefix to choose between them.

2. Sort by measure before interpolating

np.interp requires increasing x. Sorting the measure and distance arrays together handles routes digitised against the direction of numbering.

3. Interpolate the distance, then the position

Measure to distance by np.interp, distance to coordinate by LineString.interpolate. Both are piecewise linear along the vertices, which is exactly the calibration's claim.

4. Return nothing when there is nothing

A measure in a gap has no position. A nearest-segment fallback turns a data problem into a plausible wrong answer.

5. Handle from-and-to measures as a pair

A line event needs both ends on the same segment, or splitting across the segments it spans.

6. Validate against anything published

If the source also publishes coordinates, locate a sample both ways and look at the distribution of the differences.

A measure falling inside an equation gap, where no segment carries it. postmile 8.15 is real road and not a valid reference0.000 – 8.1088.198 – 19.538postmile 8.150.090 of a mile that the numbering skipsthe correct answer is None — clamping would pile every event in the gaponto the same segment boundary spatialworkflow.io
Postmile 8.15 is a real road and not a valid reference: the numbering skips it.

Placing a whole event table and checking it

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 — placing a whole event table, and checking it

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

feats = query("Traffic_AADT", "CNTY='SLO' AND RTE='101'")
aadt = pd.DataFrame([f["attributes"] for f in feats])
aadt["published"] = [Point(f["geometry"]["x"], f["geometry"]["y"]) for f in feats]
aadt["PM"] = pd.to_numeric(aadt.PM, errors="coerce")
aadt["prefix"] = aadt.PM_PFX.fillna("")
stations = aadt.drop_duplicates(subset=["prefix", "PM", "DESCRIPTION"]).dropna(subset=["PM"])

placed = [point_at(shn, pm, p) for pm, p in zip(stations.PM, stations.prefix)]
ok = [i for i, p in enumerate(placed) if p is not None]
print(f"{len(ok)} of {len(stations)} stations placed from their postmile")

mine = gpd.GeoSeries([placed[i] for i in ok], crs=shn.crs)
theirs = gpd.GeoSeries([stations.published.iloc[i] for i in ok], crs=4326).to_crs(shn.crs)
gap = mine.distance(theirs, align=False)
print(f"distance from the published point: median {gap.median():.1f} m, "
      f"90th {gap.quantile(0.9):.1f} m, max {gap.max():.1f} m")
off = theirs.distance(shn.geometry.union_all())
print(f"published points off this centreline: median {off.median():.1f} m, max {off.max():.1f} m")
48 of 48 stations placed from their postmile
distance from the published point: median 23.2 m, 90th 37.2 m, max 59.5 m
published points off this centreline: median 22.0 m, max 58.2 m

Every station is placed, and every placed point lies on the centreline; the published ones do not. They sit a median 22.0 metres to the side of it, which accounts for most of the 23.2-metre median gap, so the two disagree across the road far more than along it. Caltrans's count points used to fall on this centreline exactly; the copy saved on 9 October 2026 puts them beside the road, so this check measures where the count points are drawn, not how well a postmile places one. Test it against the layer you will actually join to.

Example 2 — a line event that spans several segments

def segment_pieces(frame, from_pm, to_pm, prefix=""):
    """Every piece of the route between two measures, segment by segment."""
    lo, hi = sorted((from_pm, to_pm))
    out = []
    for seg in frame.sort_values("bPM").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:
            out.append((round(a, 3), round(b, 3), round(b - a, 3)))
    return out

for span in [(2.0, 6.0), (7.0, 10.0), (0.0, 69.5)]:
    pieces = segment_pieces(shn, *span)
    covered = sum(p[2] for p in pieces)
    print(f"postmiles {span[0]}–{span[1]}: {len(pieces)} piece(s), "
          f"{covered:.3f} of {span[1] - span[0]:.3f} miles covered")
    for p in pieces[:3]:
        print(f"    {p[0]:7.3f} – {p[1]:7.3f}  ({p[2]:.3f} mi)")
postmiles 2.0–6.0: 1 piece(s), 4.000 of 4.000 miles covered
      2.000 –   6.000  (4.000 mi)
postmiles 7.0–10.0: 2 piece(s), 2.910 of 3.000 miles covered
      7.000 –   8.108  (1.108 mi)
      8.198 –  10.000  (1.802 mi)
postmiles 0.0–69.5: 6 piece(s), 59.429 of 69.500 miles covered
      0.000 –   8.108  (8.108 mi)
      8.198 –  19.538  (11.340 mi)
     25.825 –  32.976  (7.151 mi)

The 7.0–10.0 request comes back short: a tenth of a mile of it lies in the equation gap at 8.108, and no segment carries those measures. Reporting the covered length alongside the pieces is what makes that visible.

Example 3 — batching it

def place_many(frame, events, pm_col="PM", prefix_col="prefix"):
    index = []
    for seg in frame.itertuples():
        m = np.asarray(seg.m)
        d = along_distances(seg.geometry)
        order = np.argsort(m)
        index.append(((seg.PMPrefix or ""), m[order], d[order], seg.geometry))
    out = []
    for pm, prefix in zip(events[pm_col], events[prefix_col]):
        hit = None
        for p, m, d, geom in index:
            if p == prefix and m[0] - 1e-9 <= pm <= m[-1] + 1e-9:
                hit = geom.interpolate(float(np.interp(pm, m, d)))
                break
        out.append(hit)
    return out

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["prefix"] = ""
bridges = bridges.dropna(subset=["PM"])

pts = place_many(shn, bridges)
missing = [b for b, p in zip(bridges.itertuples(), pts) if p is None]
print(f"{len(bridges) - len(missing)} of {len(bridges)} bridges placed")
print(f"unplaced postmiles: {sorted(round(b.PM, 2) for b in missing)[:8]}")
93 of 93 bridges placed
unplaced postmiles: []

All 93 place, because every bridge on this route sits below postmile 56 and the prefixed segments start at 67.456. Run the same code on a route whose structures fall inside a realigned stretch and the blank prefix column becomes the failure: the inventory's postmiles are real, and the assumption that a missing prefix means "no prefix" is what breaks.

Why the prefix belongs in the selection

Why the prefix has to be part of the selection

Two segments in one county can cover the same postmile range under different prefixes. Selecting on the number alone picks whichever comes first in the frame, which is an arbitrary function of the download order — the worst kind of bug, because it is stable within a session and changes between them.

Why a gap must return nothing

The alternative is to clamp to the nearest valid measure, which places the event at the segment boundary. On this route that is up to a tenth of a mile from where the record says, and every event in the gap piles up at the same point, producing a cluster that looks like a real concentration of assets.

Why interpolate is the right final step

LineString.interpolate(distance) walks the vertices and returns the position at that distance, including inside a segment between two vertices. It is the exact inverse of project, so locating a point and then placing it again returns the original position — a property worth testing whenever the two functions are used together.

Distance between positions placed from postmiles and the published station coordinates. placed from the register, checked against the publishermedian23.2 m90th percentile37.2 mmaximum59.5 moff the centreline22.0 m median: mostly sideways030 m60 m48 of 48 stations placedThe published count points sit beside the road, so the gap is mostly sideways spatialworkflow.io
Placed from the register against published coordinates: a median 23.2 m, mostly sideways, because the published count points sit beside the centreline.

Boundary measures, zero-length segments and rounding

  • Measures exactly at a boundary belong to two segments; pick by prefix, then by the lower segment, and write the rule down.
  • Zero-length segments break the interpolation; filter them before building the index.
  • Rounded measures in the event table limit the precision: a hundredth of a mile is about eight metres.
  • Offsets — "30 feet right of centreline" — are a separate field, applied after the position is found.
  • PostGIS ST_LocateAlong does this against an M ordinate and returns a point with the measure attached.
  • Cache the index if you are placing many events; rebuilding the distance arrays per event dominates the runtime.

FAQ

How do I convert a milepost to a coordinate?

Find the segment whose measure range contains it, interpolate the distance along from the measures, then interpolate the position at that distance.

What if no segment contains the measure?

Return nothing. The measure is in a numbering gap, under a different prefix, or off the route.

Why does my event land in the wrong place?

Usually the prefix or the alignment was ignored when selecting the segment.

Can I use proportion instead of interpolation?

Only on an uncalibrated route; on a real one it is wrong by tens of feet.

How accurate is the result?

As accurate as the measure it came from — typically a hundredth of a mile, about eight metres — plus the centreline's own accuracy.

How do I place thousands of events?

Precompute the distance and measure arrays per segment once, then interpolate per event.