Linear Referencing Along Lines in PyQGIS

Roads, rivers, pipelines, railways and power lines are managed by distance along them, not by coordinates. A pothole is at chainage 1,240 m on route A61; a culvert sits 3.2 km downstream of a gauging station; a speed restriction runs from kilometre 12.4 to 13.1. Converting between those positions and map geometry is linear referencing, and the QGIS geometry API has the primitives for it: place a point at a distance, measure the distance of a point, and cut the piece of line between two distances.

This recipe belongs to Geometry Operations and Spatial Predicates in PyQGIS. It covers the three primitives, building an events table into geometry, generating regularly spaced points, and the checks for line direction and multipart geometry that decide whether a chainage means anything.

Distance to point, point to distance, distances to sectionA winding line from start at 0 m to end at 2,000 m. A green point sits at 600 m, placed by interpolate. An off-line orange point is projected onto the line at 1,150 m by lineLocatePoint, with a dashed perpendicular. A thick blue section runs from 1,400 m to 1,800 m, cut with a line substring.Three questions about one line0 m2,000 minterpolate(600)lineLocatePoint → 1,150substring 1,400 – 1,800

Prerequisites

  • QGIS 3.40 LTR or newer, or the QGIS 4 series.
  • A line layer in a projected CRS whose units are metres. Distances along a line in degrees are not chainages.
  • Lines that are single-part and digitised in the direction chainage increases. The checks below confirm both.

Place a point at a distance

QgsGeometry.interpolate(distance) returns a point geometry at that distance along the line, measured from the first vertex. interpolateAngle(distance) returns the line's direction there, which is useful for rotating symbols or offsetting perpendicular to the route.

import math
from qgis.core import QgsProject, QgsGeometry, QgsPointXY

routes = QgsProject.instance().mapLayersByName("road_centrelines")[0]
route = next(routes.getFeatures("\"route_id\" = 'A61'"))
line = route.geometry()

print("length:", round(line.length(), 1), "m")

chainage = 1240.0
point = line.interpolate(chainage)
angle = line.interpolateAngle(chainage)
print(point.asPoint(), "bearing", round(math.degrees(angle), 1))

offset = 4.0
normal = angle - math.pi / 2
p = point.asPoint()
kerb = QgsPointXY(p.x() + offset * math.sin(normal), p.y() + offset * math.cos(normal))

Breakdown: interpolate returns an empty geometry — not an error — when the distance is negative or longer than the line, so check point.isEmpty() when chainages come from outside data. The angle is in radians, clockwise from north, which is the convention map symbol rotation also uses. The perpendicular offset turns a chainage into a position on the left of the direction of travel: four metres to the left places an event at the kerb rather than on the centreline, which matters when a survey records which side of the road an asset is on. Subtract the right angle for the left side and add it for the right.

Measure how far along a point lies

The reverse operation projects a point onto the nearest position on the line and returns the distance from the start. Points need not lie on the line; the result is the distance of their perpendicular foot.

Measure against the right routeA main road runs left to right with a side road branching off. An asset point sits between them, slightly closer to the side road. Projecting against the nearest line gives a chainage on the side road, which is wrong. Filtering to the route the asset belongs to, then projecting, gives the correct chainage on the main road, along with the offset distance to check.Nearest line is not always the asset’s linenearest linemeasured on the side roadroute chosen by idchainage on A61, offset 18 m

assets = QgsProject.instance().mapLayersByName("drainage_assets")[0]

route_lines = {
    f["route_id"]: f.geometry() for f in routes.getFeatures()
}

results = []
for asset in assets.getFeatures():
    line = route_lines.get(asset["route_id"])
    if line is None:
        continue
    pt = asset.geometry()
    distance_along = line.lineLocatePoint(pt)
    foot = line.interpolate(distance_along)
    offset = pt.distance(foot)
    results.append((asset["asset_ref"], round(distance_along, 1), round(offset, 1)))

for ref, ch, off in sorted(results, key=lambda r: r[1]):
    flag = "  CHECK" if off > 25 else ""
    print(f"{ref:<10} ch {ch:>8} m  offset {off:>5} m{flag}")

Breakdown: Matching each asset to its own route by identifier, rather than to whichever line is nearest, is what keeps an asset at a junction on the right road; the nearest-line approach is only safe when routes never come close. lineLocatePoint returns -1 on failure, typically for an empty or non-line geometry. Interpolating back to the foot and measuring the offset is the quality check: an asset tens of metres from its route has the wrong route id or the wrong coordinates, and a chainage computed for it is fiction. Sorting by chainage gives the order a maintenance crew would drive.

Cut the section between two chainages

Linear events — a speed limit, a resurfacing scheme, a length of pipe of one material — are stored as a start and end distance. Turning them into line geometry means cutting the substring between those distances.

def section(line_geom, start, end):
    if line_geom.isMultipart():
        raise ValueError("multipart line; merge or explode first")
    start, end = sorted((max(0.0, start), min(line_geom.length(), end)))
    return QgsGeometry(line_geom.constGet().curveSubstring(start, end))

a61 = route_lines["A61"]
restriction = section(a61, 12400, 13100)
print("section length:", round(restriction.length(), 1))

Breakdown: curveSubstring is a method of the underlying curve object (QgsLineString or QgsCompoundCurve), which is why the code reaches through constGet(); wrapping the result in QgsGeometry gives back the familiar API. Clamping to the line length stops an events table with a slightly overlong end chainage — common when chainages were measured on an older version of the centreline — from producing an empty result. Sorting the two distances tolerates events recorded against the direction of the line. The multipart check matters because a substring of a MultiLineString is undefined: merge contiguous parts with lineMerge() first, or split the route.

Build geometry from an events table

The usual input is a table — a spreadsheet or database table with a route id, a start and end chainage, and attributes. The Processing algorithm does the substring for every row, but a short loop with the functions above gives you control over clamping and reporting.

From rows of chainages to line featuresA table with three rows: A61 from 12,400 to 13,100 with limit 30, A61 from 15,000 to 15,800 with limit 40, and B6154 from 200 to 26,000 with limit 20. The first two become line features on route A61. The third has an end chainage beyond the route's 4,950 m length, so it is clamped and reported for review.Every row becomes a section, or a report linespeed_limits.csvroute start end mphA61 12400 13100 30A61 15000 15800 40B6154 200 26000 202 line features30 mph40 mph1 clamped and reportedend 26,000 m beyond route length 4,950 m

from qgis.core import QgsVectorLayer, QgsFeature, QgsField, QgsFields
from qgis.PyQt.QtCore import QMetaType

events = QgsProject.instance().mapLayersByName("speed_limits")[0]   # geometryless table

out = QgsVectorLayer(f"LineString?crs={routes.crs().authid()}", "speed_limit_sections", "memory")
out_fields = QgsFields()
for name, kind in (("route_id", QMetaType.Type.QString), ("start_m", QMetaType.Type.Double),
                   ("end_m", QMetaType.Type.Double), ("mph", QMetaType.Type.Int)):
    out_fields.append(QgsField(name, kind))
out.dataProvider().addAttributes(out_fields)
out.updateFields()

features, report = [], []
for ev in events.getFeatures():
    line = route_lines.get(ev["route_id"])
    if line is None:
        report.append(f"{ev['route_id']}: no such route")
        continue
    if ev["end_m"] > line.length() + 1:
        report.append(f"{ev['route_id']} {ev['start_m']}–{ev['end_m']}: "
                      f"beyond length {line.length():.0f}")
    f = QgsFeature(out_fields)
    f.setGeometry(section(line, ev["start_m"], ev["end_m"]))
    f.setAttributes([ev["route_id"], ev["start_m"], ev["end_m"], ev["mph"]])
    features.append(f)

out.dataProvider().addFeatures(features)
QgsProject.instance().addMapLayer(out)
print("\n".join(report) or "all events within route lengths")

Breakdown: The one-metre tolerance on the length check avoids reporting rounding differences, while still catching events recorded against a different, longer version of the route. Clamped events are still created — the report is for review, not rejection — because a partly wrong section is usually more useful on the map than a missing one. The memory layer is the quick output; save it with native:savefeatures when the result is final. For a table that is already clean, native:linesubstring in Processing, or a virtual layer as in querying layers with virtual layer SQL, are alternatives.

Direction, multipart lines and evenly spaced points

A chainage is meaningless if the line runs the wrong way. Before trusting any result, confirm the first vertex is where chainage zero should be, and reverse lines that are backwards.

start_node = QgsGeometry.fromPointXY(QgsPointXY(430112.0, 431050.0))  # known 0 m marker

for rid, geom in route_lines.items():
    first = QgsGeometry(geom.vertexAt(0))
    if geom.isMultipart():
        print(rid, "is multipart:", geom.constGet().numGeometries(), "parts")
    if rid == "A61" and first.distance(start_node) > 50:
        route_lines[rid] = QgsGeometry(geom.constGet().reversed())
        print(rid, "reversed so chainage starts at the known marker")

markers = processing.run("native:pointsalonglines", {
    "INPUT": routes, "DISTANCE": 100, "START_OFFSET": 0, "END_OFFSET": 0,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]
print(markers.featureCount(), "100 m markers; 'distance' field holds each chainage")

Breakdown: Comparing the first vertex against a known zero marker is the only reliable direction test, because direction is a convention of the organisation that owns the route, not a property of the geometry. reversed() flips a curve's vertex order without changing its shape. native:pointsalonglines generates a point every 100 m and writes the chainage of each into a distance field, which is a ready-made set of marker posts for labelling, and the way to sample a raster along a route as in extracting an elevation profile. Remember to import processing first when running outside the Python console.

QGIS version compatibility

interpolate, interpolateAngle and lineLocatePoint have been part of QgsGeometry since QGIS 3.0; curveSubstring since 3.4. QgsField with QMetaType.Type is the 3.38+ constructor and the only one on the QGIS 4 series; on older releases use QVariant.String, QVariant.Double and QVariant.Int. Expression equivalents — line_interpolate_point, line_locate_point, line_substring — are available from 3.4 for use in field calculations and virtual fields.

Troubleshooting

  • interpolate returns an empty geometry. The distance is beyond the line length or negative.
  • Chainages run backwards. The line is digitised end to start; reverse it.
  • Results on the wrong road near junctions. Points were measured against the nearest line instead of their own route.
  • Substrings are empty for multipart routes. Merge parts with lineMerge() or explode and measure part by part with cumulative offsets.
  • Distances look like tiny decimals. The layer is in degrees; reproject to a metric CRS first.

Conclusion

Keep routes single-part, projected in metres and digitised in chainage order, and check all three before trusting a result. Use interpolate to place events, lineLocatePoint with the right route and an offset check to measure them, and curveSubstring to cut sections from start and end distances — clamping and reporting rather than silently dropping the rows that do not fit.

Frequently Asked Questions

Can I use M values instead of computing distances? Yes, if the lines carry calibrated measures. Read them per vertex with QgsPoint.m() and interpolate between vertices; this is the right approach when official chainages do not match geometric length.

How do I handle routes that are longer on the ground than on the map? Calibrate: scale geometric distance by the ratio of official length to geometric length, or store M values from known marker posts.

Is there a Processing algorithm for events tables?native:linesubstring cuts one distance range for every feature of a line layer; combine it with a join for per-row ranges, or use the loop above.

Does this work on curved geometries? Yes. curveSubstring and interpolate work on compound curves and circular strings.