How to read a linear referencing system in Python

On this page

Problem statement

Route data with measures arrives from a feature service, and three things go wrong before any analysis starts: the service returns 2,000 features and stops without saying so, the measures are absent unless you ask for them, and a multipart feature's attributes describe all of its parts at once.

This guide downloads the California State Highway Network with its postmiles intact, turns it into a GeoDataFrame where the measures survive reprojection, and stores it in a format that keeps them.

Quick answer

Ask for M, page past the record limit, and keep one row per path:

import hashlib, urllib.parse
import geopandas as gpd, numpy as np, requests
from shapely.geometry import LineString

BASE = ("https://caltrans-gis.dot.ca.gov/arcgis/rest/services/CHhighway/"
        "SHN_Lines/FeatureServer/0/query")
# each page as Caltrans served it on 9 October 2026, named by a SHA-1 of the live URL
FROZEN = "https://data.spatialworkflow.io/caltrans-2026-10-09/{key}.json"

def fetch(where, out_sr=4326):
    features, offset = [], 0
    while True:
        params = {"where": where, "outFields": "*", "f": "json", "outSR": out_sr,
                  "returnM": "true",                      # without this there is no measure
                  "resultOffset": offset, "resultRecordCount": 2000}
        url = BASE + "?" + 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 tells you
            return features
        offset += len(page["features"])

feats = fetch("District=5")
paths = sum(len(f["geometry"]["paths"]) for f in feats)
print(f"{len(feats)} features, {paths} paths, "
      f"{sum(len(p) for f in feats for p in f['geometry']['paths']):,} vertices")
print(f"multipart features: {sum(1 for f in feats if len(f['geometry']['paths']) > 1)}")
print(f"a vertex looks like: {np.round(feats[0]['geometry']['paths'][0][0], 4).tolist()}")

480 features, 482 paths, 63,216 vertices
multipart features: 2
a vertex looks like: [-121.6665, 36.6385, 18.737]

Three numbers per vertex: X, Y and the postmile. Without returnM the third is simply missing, and nothing in the response says the measures were dropped.

Map of Caltrans District 5 on the central California coast with every state highway coloured by its postmile, from yellow near zero to dark purple above 60; the colours restart at yellow at each county boundary. — The third number of every vertex, drawn: District 5’s state highways coloured by postmile. Each road starts again at zero — yellow — where it enters a county, which is why US-101 changes from purple back to yellow at each county line. A postmile names a place only together with its county and route.

Four things to get right between the service and the dataframe.

Filter on the server, request measures, page to the end

1. Filter on the server

where clauses run in the database. Downloading a state to select a county is slow and unnecessary.

2. Ask for the measures explicitly

returnM=true for Esri services; hasM on the layer description says whether there are any.

3. Page until the service says it is done

exceededTransferLimit is the flag. A loop that stops at the first short page is right; a loop that assumes one request is enough is silently wrong at exactly 2,000 features.

4. Expand multipart features into rows

One geometry per row makes every later operation simpler, and the measures belong to the path rather than to the feature.

5. Separate the measures from the geometry

Keep the measures in their own column. They survive to_crs, explode, joins and Parquet; an M ordinate does not survive the first of those.

6. Store it where the measures come back

GeoParquet keeps a measures column, and keeps an M ordinate if you have one. Most OGR formats do neither.

The 2,000-record limit is not an error; it is a flag you have to read.

A GeoDataFrame whose measures survive

Example 1 — a GeoDataFrame with measures that survive

def to_frame(features, out_sr=4326):
    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])
    g = gpd.GeoDataFrame(rows, geometry=geoms, crs=out_sr)
    g["m"] = measures
    return g

shn = to_frame(feats)
shn["key"] = [f"{r.County}-{r.Route}{r.PMPrefix or ''}{r.PMSuffix or ''}-{r.AlignCode}"
              for r in shn.itertuples()]

metric = shn.to_crs(3310)                     # measures unaffected: they are a column
print(f"{len(metric)} rows, {metric.key.nunique()} distinct route keys")
print(f"total centreline: {metric.geometry.length.sum() / 1609.344:,.0f} miles")
print(metric.groupby("AlignCode").size().to_string())
print(f"measures still here after reprojection: "
      f"{metric.m.iloc[0][:3].round(3).tolist()} …")

482 rows, 186 distinct route keys
total centreline: 2,298 miles
AlignCode
Left                 236
Left Independent       4
Right                238
Right Independent      4
measures still here after reprojection: [18.737, 18.843, 18.925] …

Example 2 — checking what came back

import pandas as pd

checks = pd.DataFrame({
    "rows": [len(shn)],
    "with measures": [(shn.m.map(len) > 1).sum()],
    "monotonic": [sum(bool(np.all(np.diff(m) >= 0) or np.all(np.diff(m) <= 0))
                      for m in shn.m)],
    "increasing": [sum(bool(np.all(np.diff(m) >= 0)) for m in shn.m)],
    "path span == attribute span": [
        sum(abs(abs(m[-1] - m[0]) - abs(r.ePM - r.bPM)) < 0.002
            for m, r in zip(shn.m, shn.itertuples()))],
})
print(checks.T.rename(columns={0: "count"}).to_string())

backwards = [r for m, r in zip(shn.m, shn.itertuples()) if m[0] > m[-1]]
print(f"\npaths whose measures decrease along the geometry: {len(backwards)}")
if backwards:
    r = backwards[0]
    print(f"  e.g. {r.County}-{r.Route} {r.AlignCode}: "
          f"M runs {r.m[0]:.3f} -> {r.m[-1]:.3f}")

                             count
rows                           482
with measures                  482
monotonic                      482
increasing                     482
path span == attribute span    478

paths whose measures decrease along the geometry: 0

Every path in this district increases, and four have a measure span that does not match their feature's bPM–ePM: those are the two multipart features. Decreasing measures are legal — they mean the geometry was digitised against the direction of increasing postmile — so the check is worth running rather than assuming, and every interpolation has to handle both directions.

Example 3 — storing it

out = metric[["County", "Route", "PMPrefix", "AlignCode", "bPM", "ePM",
              "bOdometer", "eOdometer", "key", "m", "geometry"]]
out.to_parquet("/tmp/shn_d5.parquet")

back = gpd.read_parquet("/tmp/shn_d5.parquet")
print(f"round trip: {len(back)} rows, crs {back.crs.to_epsg()}, "
      f"measures type {type(back.m.iloc[0]).__name__} of length {len(back.m.iloc[0])}")
print(f"first three measures: {np.asarray(back.m.iloc[0])[:3].round(3).tolist()}")

out.to_file("/tmp/shn_d5.gpkg", driver="GPKG")
gpkg = gpd.read_file("/tmp/shn_d5.gpkg")
print(f"GeoPackage kept a column called m: {'m' in gpkg.columns}, "
      f"as a {type(gpkg.m.iloc[0]).__name__}")
print(f"  {str(gpkg.m.iloc[0])[:60]}…")

round trip: 482 rows, crs 3310, measures type ndarray of length 97
first three measures: [18.737, 18.843, 18.925]
GeoPackage kept a column called m: True, as a str
  [18.737  18.8427 18.9248 19.1387 19.1864 19.19   19.3549 19.…

GeoPackage does not refuse the array column — it writes str() of it, complete with numpy's abbreviating ellipsis, so the file contains a truncated printout of the measures that will parse back as text. That is worse than an error, and it is why Parquet is the right store here: OGR has no column type for an array, so the alternatives are one row per vertex, an explicit JSON string, or silent damage.

Why the transfer-limit flag matters

Why `exceededTransferLimit` matters more than it looks

The response with 2,000 features is a valid, complete-looking GeoJSON. Nothing about it says "there is more". Any script that takes the first page and moves on produces an analysis of an arbitrary subset, and the subset changes whenever the service's ordering does.

Why one row per path

A multipart feature's attributes apply to the feature. Its measures apply to each path separately, so bPM and ePM bound the union of the parts and match none of them. Exploding first makes the row the unit of both geometry and measure, which is the only arrangement where the arithmetic is simple.

Why the measures go in a column

An M ordinate is the standards-correct place for them and it does not survive a reprojection in this stack. A column of arrays survives everything and costs one line of bookkeeping — keep them the same length as the coordinate array and every operation stays obvious.

Where the measures can live, and what each choice costs.

Server-side reprojection, precision and rate limits

outSR on the request reprojects server-side and keeps M; doing it locally does not.
Geometry precision: request geometryPrecision to cut payload size where the vertices are over-precise.
Layer metadata at .../FeatureServer/0?f=json gives hasM, maxRecordCount and the field list.
Rate limits: page with a small delay on public services, and cache the result.
Curved segments in an Esri service can come back as true curves; f=json with returnTrueCurves off densifies them.
OBJECTID is not stable between versions of the layer; key on the route fields.

Internal links

Linear referencing explained: routes, measures and events — the model behind the fields
M values explained: geometry that carries a measure — why the column, not the ordinate
How to find the measure of a point on a route — the first thing to do with it
How to paginate a large feature API in Python — paging in general
GeoJSON downloaded from an API is empty or truncated — the failure this avoids
Route calibration explained — what the measures mean

FAQ

How do I get M values from an Esri feature service?

Pass returnM=true; the measures arrive as a third number on each vertex.

How do I know the download is complete?

Check exceededTransferLimit in the response and page with resultOffset until it is false.

Should I explode multipart routes?

Yes. Measures belong to the path, and the feature's begin and end measures cover all the parts together.

Why keep measures in a column rather than on the geometry?

Because reprojection drops the M ordinate, and most formats will not write it.

Which format should I store routes in?

GeoParquet — it keeps both the geometry and an array column of measures.

Can I filter by postmile on the server?

Yes, with a where clause on the begin and end postmile fields, but mind the prefix and the county.