How to find the measure of a point on a route

Problem statement

Something has a coordinate and needs a measure: a GPS-logged sign, a crash imported from a mapping app, an asset digitised from imagery. The operation is "project the point onto the route and read the measure at that position", and each of its three parts has a way to go wrong — the wrong route, the wrong projection, and a measure read from proportion rather than from the calibration.

This guide locates points on the California State Highway Network, checks the answers against published postmiles, and shows how far off the two carriageways of a divided highway put you.

Quick answer

Project in a metric CRS, then interpolate the measure at that distance:

import numpy as np
from shapely.geometry import Point
from lrs import routes

shn = routes("County='SLO' AND Route=101")        # both alignments, metric CRS

def measure_at(line, m, point):
    """Interpolate against distance along the line, which always increases,
    so measures that run backwards are handled without a special case."""
    xy = np.asarray(line.coords)
    along = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))])
    return float(np.interp(line.project(point), along, m))

here = Point(*np.asarray(shn.geometry.iloc[0].coords)[len(shn.m.iloc[0]) // 2])
for seg in shn.itertuples():
    gap = seg.geometry.distance(here)
    if gap < 50:
        print(f"{seg.County}-{seg.Route}{seg.PMPrefix or ''} {seg.AlignCode:6} "
              f"{gap:6.1f} m away, postmile {measure_at(seg.geometry, seg.m, here):.3f}")

SLO-101R Right     0.0 m away, postmile 68.483
SLO-101R Left     17.0 m away, postmile 68.468

The point is within fifty metres of both carriageways, and each returns its own postmile. Choosing between them is part of the operation, not an afterthought.

Projected CRS, narrowed candidates, project then interpolate

1. Work in a projected CRS

project and distance return numbers in whatever units the coordinates are in. In degrees those numbers are meaningless.

2. Narrow the candidates before projecting

A spatial index, a buffer, a route number from the source data — anything that stops you projecting onto 482 segments and taking the nearest.

3. Project, then interpolate

line.project(point) gives the distance along; np.interp against the stored measures turns that into a measure. Skipping the second step gives you a proportion, not a postmile.

4. Keep the offset

The distance from the point to the line is the quality indicator. Twenty metres is a carriageway; two hundred is a different road.

5. Resolve the alignment deliberately

Direction of travel, a route field on the source, or the side the point falls on. Not "whichever is nearer" — on a divided highway that flips back and forth along the route.

6. Report the unmatched

Points with no candidate within a sensible distance are not on the network, and should come back as null rather than as the nearest thing found.

Locating and scoring a real inventory

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 — locating a real inventory and scoring it

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

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 = bridges.dropna(subset=["PM", "LAT", "LON"])
pts = gpd.GeoSeries([Point(x, y) for x, y in zip(bridges.LON, bridges.LAT)],
                    crs=4326).to_crs(shn.crs)

rows = []
for (_, b), p in zip(bridges.iterrows(), pts):
    best = min(((seg.geometry.distance(p), seg) for seg in shn.itertuples()),
               key=lambda t: t[0])
    offset, seg = best
    rows.append({"bridge": b.BRIDGE, "published_pm": b.PM,
                 "located_pm": measure_at(seg.geometry, seg.m, p),
                 "align": seg.AlignCode, "offset_m": offset})
located = pd.DataFrame(rows)
located["error_pm"] = (located.located_pm - located.published_pm).abs()

print(f"{len(located)} bridges located")
print(f"offset to the nearest carriageway: median {located.offset_m.median():.1f} m, "
      f"90th {located.offset_m.quantile(0.9):.1f} m")
print(f"postmile error against the published value: median "
      f"{located.error_pm.median():.3f} mi, 90th {located.error_pm.quantile(0.9):.3f} mi, "
      f"max {located.error_pm.max():.3f} mi")
print(located.nlargest(4, "error_pm")[["bridge", "published_pm", "located_pm",
                                       "align", "offset_m"]].to_string(index=False))

93 bridges located
offset to the nearest carriageway: median 5.9 m, 90th 19.1 m
postmile error against the published value: median 0.006 mi, 90th 0.023 mi, max 0.087 mi
  bridge  published_pm  located_pm align  offset_m
49 0060L         35.60   35.513325  Left  5.236057
 49 0007         36.58   36.532004 Right  6.460902
 49 0224         55.97   55.927192  Left 77.961141
49 0252L          6.45    6.411590  Left  2.270114

A median error of 0.006 miles is 32 feet — the published postmile rounded to a hundredth of a mile, plus the structure's offset from the centreline. The worst case is a bridge 78 metres from the nearest carriageway: an overcrossing, whose coordinate is on the road passing above rather than on US 101, and whose large offset is exactly the signal that says so.

Example 2 — what the alignment choice costs

per_align = {}
for align in ("Left", "Right"):
    sub = shn[shn.AlignCode == align]
    out = []
    for (_, b), p in zip(bridges.iterrows(), pts):
        offset, seg = min(((s.geometry.distance(p), s) for s in sub.itertuples()),
                          key=lambda t: t[0])
        out.append(abs(measure_at(seg.geometry, seg.m, p) - b.PM))
    per_align[align] = np.array(out)
    print(f"forced onto the {align:5} carriageway: median error "
          f"{np.median(per_align[align]):.3f} mi, max {per_align[align].max():.3f} mi")

flip = sum(1 for (_, b), p in zip(bridges.iterrows(), pts)
           if (min(shn[shn.AlignCode == 'Left'].geometry.distance(p)) <
               min(shn[shn.AlignCode == 'Right'].geometry.distance(p))))
print(f"bridges whose nearest carriageway is the Left one: {flip} of {len(bridges)}")

forced onto the Left  carriageway: median error 0.006 mi, max 0.087 mi
forced onto the Right carriageway: median error 0.006 mi, max 0.088 mi
bridges whose nearest carriageway is the Left one: 50 of 93

The two carriageways give almost the same postmile on this route, because Caltrans numbers them together. On a system that numbers each direction independently the same experiment would show a systematic offset, and "nearest line wins" would scatter events between the two.

Example 3 — a fast lookup for many points

from shapely import STRtree

tree = STRtree(shn.geometry.values)
sample = pts.sample(min(len(pts), 60), random_state=0)

def locate(point, max_offset=100.0):
    idx = tree.query_nearest(point, max_distance=max_offset)
    if len(idx) == 0:
        return None
    seg = shn.iloc[int(idx[0])]
    return {"key": f"{seg.County}-{seg.Route}{seg.PMPrefix or ''}-{seg.AlignCode}",
            "pm": measure_at(seg.geometry, seg.m, point),
            "offset_m": seg.geometry.distance(point)}

found = [locate(p) for p in sample]
hits = [f for f in found if f]
print(f"{len(hits)} of {len(sample)} points within 100 m of a route")
print(f"median offset {np.median([h['offset_m'] for h in hits]):.1f} m")
for h in hits[:4]:
    print(f"  {h['key']:22} postmile {h['pm']:7.3f}  offset {h['offset_m']:5.1f} m")

60 of 60 points within 100 m of a route
median offset 5.8 m
  SLO-101-Left           postmile  49.394  offset   3.6 m
  SLO-101-Left           postmile   0.801  offset   5.6 m
  SLO-101-Right          postmile  65.056  offset  19.3 m
  SLO-101-Right          postmile  28.808  offset   7.8 m

STRtree.query_nearest with a maximum distance does the candidate narrowing and the "not on the network" test in one call, which is what makes this scale past a few hundred points.

Why the offset is part of the answer

Why the offset is part of the answer

A projection always succeeds. The only evidence that the result means anything is how far the point was from the line, so an LRS location without its offset cannot be audited. Store it, and set a threshold that reflects the network: a hundred metres on rural highway, ten in a city where parallel streets are close.

Why "nearest" is not a rule

Nearest is a sensible default and a bad rule. Two carriageways, a frontage road, an on-ramp and the mainline are all within a few tens of metres of each other at an interchange, and the nearest one changes between adjacent points. Where the source data names a route, use it; where it names a direction, use that to pick the alignment.

Why the measure must come from interpolation

project returns a distance along the geometry. The measure at that distance is whatever the calibration says, and on a real route the two are not proportional — a point two thirds of the way along a segment is not at two thirds of its postmile range.

Points past the ends and overlapping routes

• Points beyond the ends project to the endpoint, and project gives 0 or the full length — check before trusting the measure.
• Self-overlapping routes (a loop, a switchback) have two positions at the same distance from a point; projection picks one.
• Very long segments make project slow; split them or index them.
• A point exactly between two carriageways is a coin flip; use the direction field.
• PostGIS has ST_LineLocatePoint and ST_InterpolatePoint, which do this against the M ordinate directly.
• Record the route version: a located event is only valid against the centreline it was located on.

Internal links

• How to read a linear referencing system in Python — getting the routes and measures
• How to find the point at a given measure on a route — the inverse operation
• A point snaps to the wrong carriageway — when the alignment goes wrong
• Route calibration explained — why the measure is not a proportion
• How to snap points to a street network in Python — the same projection for routing
• Linear referencing explained: routes, measures and events — the model

FAQ

How do I get a postmile from a coordinate?

Project the point onto the route to get a distance along, then interpolate the stored measures at that distance.

Which CRS should I use?

A projected one in metres. Projection distances in degrees are not distances.

How do I pick between two nearby routes?

Use whatever the source says — route number, direction, jurisdiction — and fall back to nearest only with the offset recorded.

What offset is acceptable?

It depends on the network: tens of metres on a highway, a few metres in a dense street network. Choose it explicitly.

Why is my measure slightly off the published one?

Published postmiles are rounded, and the structure is not on the centreline. A median error of a few hundredths of a mile is normal.

How do I locate thousands of points?

Build an STRtree over the routes and use query_nearest with a maximum distance.