🔧 Fix Linear referencing Files & formats Linear Referencing · step 19 Notebook By · Updated

M values are lost when you export the route

On this page

Problem statement

The routes came out of the service with a postmile on every vertex. After a reprojection and a write to GeoPackage, the measures are gone — or worse, they are in a column that looks fine and holds a printout of a numpy array, which on a long route is truncated with an ellipsis.

Nothing raised. to_crs returned geometry, the file wrote, the read-back has the same number of rows. The measures simply are not there, and the next step in the pipeline interpolates positions from proportion instead, which produces answers that are wrong by tens of feet and look right.

Quick answer

Test what your stack actually keeps, before you build on it:

import geopandas as gpd, numpy as np, shapely
from shapely import from_wkt, has_m

route = from_wkt("LINESTRING M (-120.9 35.4 10, -120.8 35.5 20, -120.7 35.7 34)")
g = gpd.GeoDataFrame({"id": [1]}, geometry=[route], crs=4326)
print(f"shapely {shapely.__version__} / GEOS {'.'.join(map(str, shapely.geos_version))}")
print(f"as built:            M {bool(has_m(g.geometry.values)[0])}")
print(f"after to_crs(3310):  M {bool(has_m(g.to_crs(3310).geometry.values)[0])}")

g.to_parquet("/tmp/lr_keep.parquet")
print(f"GeoParquet:          M "
      f"{bool(has_m(gpd.read_parquet('/tmp/lr_keep.parquet').geometry.values)[0])}")

for driver, path in [("GeoJSON", "/tmp/m.geojson"), ("GPKG", "/tmp/m.gpkg"),
                     ("ESRI Shapefile", "/tmp/m.shp"), ("FlatGeobuf", "/tmp/m.fgb")]:
    try:
        g.to_file(path, driver=driver)
        print(f"{driver:15} wrote, M "
              f"{bool(has_m(gpd.read_file(path).geometry.values)[0])}")
    except Exception as e:
        print(f"{driver:15} {type(e).__name__}: {str(e)[:52]}")
shapely 2.1.2 / GEOS 3.13.1
as built:            M True
after to_crs(3310):  M False
GeoParquet:          M True
GeoJSON         GeometryError: Could not create geometry at index 0 for WKB type 10
GPKG            GeometryError: Could not create geometry at index 0 for WKB type 10
ESRI Shapefile  GeometryError: Could not create geometry at index 0 for WKB type 10
FlatGeobuf      GeometryError: Could not create geometry at index 0 for WKB type 10

Two failure modes in one run: a reprojection that succeeds and drops the measures, and four drivers that refuse outright. The refusal is the good case.

Operations and formats with whether each keeps, drops or refuses an M ordinate. keeps, drops, refusesoperationoutcomehow you find outto_crs()drops, silentlyhas_m() afterwardsGeoParquet round tripkeepsnothing to find outGeoJSON / GPKG / SHP / FGBrefusesGeometryError on writesimplify(), segmentize()breaks the pairinglen(m) != len(coords)WKT / WKBkeepsthe only way to construct Mthe middle column is where a pipeline loses data without an errorshapely 2.1.2, GEOS 3.13.1, GeoPandas 1.1 spatialworkflow.io
Three outcomes per operation, and the middle column is the dangerous one.

Decide where measures live, keep a column, store GeoParquet

1. Decide where the measures live, once

Either on the geometry as M, or in a column beside it. Mixing the two is how a pipeline ends up with a geometry that has lost them and a column that was never updated.

2. Prefer a measures column for anything computed

An array column survives to_crs, explode, joins, concatenation and Parquet. It is the only representation in this stack that survives all of them.

3. Use GeoParquet for storage

It round-trips both a measures column and an M ordinate, because it stores WKB and does not go through an OGR driver.

4. Convert to M only at the boundary

If a consumer wants LINESTRING M, build the WKT at export time. That is the one construction route shapely offers.

5. Assert, do not assume

One line — assert bool(has_m(...)) or assert len(m) == len(coords) — at each pipeline stage catches the silent loss where it happens rather than three steps later.

6. Check the round trip, not the write

A driver that writes without complaint may still have discarded the measures. Read the file back and test it.

An array column written to GeoPackage coming back as a printout with an ellipsis. the worst outcome is not an errornp.ndarray, 1,500 values[0.0, 1.0, 2.0, … , 1499.0]to_filestr, 65 characters[0.e+00 1.e+00 ... 1.499e+03]GeoPackage has no array column type, so it writes str() of the arraybelow numpy’s 1,000-item threshold it is complete but must be re-parsed;above it, 1,494 of the 1,500 measures are replaced by an ellipsis spatialworkflow.io
The worst outcome: a column that survives as a printout of itself, ellipsis and all.

A safe export that keeps everything

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 — the measures column, and what writes it

from lrs import routes

shn = routes("County='SLO' AND Route=101 AND AlignCode='Right'").head(3)
print(f"{len(shn)} segments, measures length matches vertices: "
      f"{[len(m) == len(g.coords) for m, g in zip(shn.m, shn.geometry)]}")

shn.to_parquet("/tmp/lr_col.parquet")
back = gpd.read_parquet("/tmp/lr_col.parquet")
print(f"Parquet round trip: type {type(back.m.iloc[0]).__name__}, "
      f"length {len(back.m.iloc[0])}, first {np.asarray(back.m.iloc[0])[:3].round(3)}")

shn.to_file("/tmp/lr_col.gpkg", driver="GPKG")
gpkg = gpd.read_file("/tmp/lr_col.gpkg")
longest = max(range(len(gpkg)), key=lambda i: len(shn.m.iloc[i]))
text = str(gpkg.m.iloc[longest])
print(f"GeoPackage round trip: type {type(gpkg.m.iloc[0]).__name__}")
print(f"  the row with {len(shn.m.iloc[longest])} vertices came back as "
      f"{len(text)} characters of whitespace-delimited text")
print(f"  {text[:52]}… {text[-20:]}")
print(f"  abbreviated by numpy: {'...' in text}   "
      f"(numpy prints arrays in full below {np.get_printoptions()['threshold']} items)")

big = gpd.GeoDataFrame({"m": [np.arange(1500.0)]}, geometry=shn.geometry.iloc[:1].values,
                       crs=shn.crs)
big.to_file("/tmp/lr_big.gpkg", driver="GPKG")
back_big = str(gpd.read_file("/tmp/lr_big.gpkg").m.iloc[0])
print(f"  a 1,500-value array: {len(back_big)} characters, "
      f"abbreviated: {'...' in back_big}")
print(f"  {back_big[:44]} … {back_big[-30:]}")
3 segments, measures length matches vertices: [True, True, True]
Parquet round trip: type ndarray, length 23, first [67.456 67.478 67.712]
GeoPackage round trip: type str
  the row with 358 vertices came back as 2894 characters of whitespace-delimited text
  [49.477  49.5357 49.635  49.6695 49.6935 49.7127 49.… 247 67.3609 67.456 ]
  abbreviated by numpy: False   (numpy prints arrays in full below 1000 items)
  a 1,500-value array: 65 characters, abbreviated: True
  [0.000e+00 1.000e+00 2.000e+00 ... 1.497e+03 … 1.497e+03 1.498e+03 1.499e+03]

The stringified column is not always truncated — numpy prints arrays in full up to a thousand items, so a 358-vertex route survives as parseable whitespace-delimited text. Above that threshold it does not: a 1,500-value array comes back with its middle replaced by a literal ..., so the file holds the first three values, the last three, and an ellipsis where the other 1,494 were. Either way the column is now text that has to be re-parsed, and whether it is complete depends on a numpy print option.

Example 2 — a safe export that keeps everything

def to_wkt_m(geom, measures):
    """LINESTRING M WKT — the one construction path shapely offers."""
    coords = np.asarray(geom.coords)
    body = ", ".join(f"{x} {y} {m}" for (x, y), m in zip(coords, measures))
    return f"LINESTRING M ({body})"

shn2 = shn.copy()
shn2["wkt_m"] = [to_wkt_m(g, m) for g, m in zip(shn2.geometry, shn2.m)]
shn2["m_json"] = [list(np.round(np.asarray(m), 4)) for m in shn2.m]
export = shn2.drop(columns=["m"])
export["m_json"] = export.m_json.map(lambda v: ",".join(map(str, v)))
export.to_file("/tmp/lr_safe.gpkg", driver="GPKG")

back = gpd.read_file("/tmp/lr_safe.gpkg")
restored = [np.array([float(x) for x in s.split(",")]) for s in back.m_json]
print(f"restored measures: {[len(r) for r in restored]} values per row")
print(f"identical to the originals: "
      f"{all(np.allclose(a, np.round(np.asarray(b), 4)) for a, b in zip(restored, shn.m))}")
print(f"the WKT column also survives: "
      f"{bool(has_m(from_wkt(back.wkt_m.iloc[0])))}, "
      f"{len(from_wkt(back.wkt_m.iloc[0]).coords)} vertices")
restored measures: [23, 358, 23] values per row
identical to the originals: True
the WKT column also survives: True, 23 vertices

Two belt-and-braces columns: a delimited string of measures that any reader can parse, and the OGC-standard WKT for anything that understands M. Both are text, both survive every driver, and neither is silently truncated.

Example 3 — the assertion that would have caught it

def with_measures(frame):
    """Every row's measures must match its vertices; raise if not."""
    bad = [i for i, (m, g) in enumerate(zip(frame.m, frame.geometry))
           if len(m) != len(g.coords)]
    if bad:
        raise ValueError(f"{len(bad)} rows have measures that do not match the geometry")
    return frame

ok = with_measures(shn)
print(f"check passed on {len(ok)} rows")

broken = shn.copy()
broken["geometry"] = broken.geometry.simplify(50)      # drops vertices, keeps m
try:
    with_measures(broken)
except ValueError as e:
    print(f"after simplify(50): {e}")
print(f"  vertices before/after: "
      f"{[len(g.coords) for g in shn.geometry]} -> "
      f"{[len(g.coords) for g in broken.geometry]}")
check passed on 3 rows
after simplify(50): 3 rows have measures that do not match the geometry
  vertices before/after: [23, 358, 23] -> [3, 28, 3]

simplify is the other silent measure-killer: it rewrites the geometry with fewer vertices and leaves the measures column at its old length, so the arrays no longer correspond. Nothing raises until something interpolates.

Why reprojection drops M

Why reprojection drops M

Transforming coordinates rebuilds the geometry from the transformed X and Y. The transformer has no concept of a measure and no reason to carry one, so the rebuilt geometry is XY. It is a defensible implementation and a silent data loss, which is why the assertion matters more than the argument.

Why a stringified array is worse than an error

GeoPackage has no array column type, so writing one falls back to str(). What lands in the file is a printout: whitespace-delimited, no declared precision, and above numpy's thousand-item threshold, abbreviated with an ellipsis. It reads back as a plausible text column, and code that parses it gets either a lossy round trip or, on a long route, a shorter array than it started with.

Why the WKT column is worth the bytes

LINESTRING M in a text column is the standard, portable, self-describing representation of measured geometry. It costs perhaps 30% more bytes than the WKB, it survives every format, and it can be read by anything that implements the OGC standard — including PostGIS, where it becomes a real measured geometry on import.

The assertions that catch measure loss where it happens. the assertions that catch ithas_m() after every reprojectionlen(measures) == len(geom.coords) per rowafter every geometry editread the file back and test ita clean write is not a kept measureno array columns written through OGRexport a LINESTRING M WKT column at the boundarytrust that the write succeededEach one catches the loss at the step that caused it, not three steps later spatialworkflow.io
Six checks, each one line, each catching a loss where it happens.

Simplify, explode, GDAL and CSV

  • simplify, segmentize, snap all change the vertex count and invalidate a parallel measures column.
  • explode keeps the column per row, so a multipart row's measures apply to all its parts — split them deliberately.
  • PostGIS stores M natively; ST_Force4D, ST_AddMeasure and ST_LocateAlong all work on it.
  • GDAL can write M to some formats from the command line even where this Python stack refuses; check your versions rather than assuming either way.
  • CSV exports of measured routes want one row per vertex, with route, sequence, x, y and m.
  • Compression: a measures column of float64 compresses well in Parquet; there is no reason to round it for space.

FAQ

Why did my M values disappear after reprojecting?

to_crs rebuilds the geometry from X and Y only. Keep measures in a column, or re-attach them afterwards.

Which formats keep M?

GeoParquet and WKB/WKT in this stack; the OGR drivers here refuse to write it at all.

What happens if I write an array column to GeoPackage?

It is written as str() of the array — whitespace-delimited text, and above numpy's thousand-item print threshold, truncated with an ellipsis.

How do I export measured routes safely?

A delimited measures column plus a LINESTRING M WKT column, both of which are text.

What breaks a measures column?

Anything that changes the vertex count: simplify, segmentize, snapping, or a manual geometry edit.

How do I catch this early?

Assert that every row's measures array has the same length as its coordinate array, after every geometry operation.