Calculating Road Curvature with Python Shapely

Compute per-vertex road curvature ($\kappa$) over discrete Shapely LineString centerlines, the step that turns vectorized lane geometry into a curvature attribute consumed by trajectory planning, superelevation modeling, and lane schema population inside a production Lane Geometry Extraction & Road Network Processing pipeline. Shapely exposes GEOS-backed geometric predicates but no differential-geometry operators, so $\kappa$ must be derived numerically from the coordinate sequence — vectorized through NumPy, projection-aware, and held to a ±0.001 m⁻¹ accuracy budget against closed-form benchmarks under a fixed per-worker RAM ceiling.

Vectorized curvature from a discrete LineString, with stability guards and a validation gate:

Curvature at any point along a planar curve is defined as:

$$\kappa = \frac{|x'y'' - y'x''|}{(x'^2 + y'^2)^{3/2}}$$

where primes denote derivatives with respect to arc length $s$. Discrete LineString geometries have no analytical derivatives, so $x'$, $y'$, $x''$, $y''$ are recovered by finite differences against a cumulative arc-length parameter.

Prerequisites #

• Python 3.10+, Shapely 2.0+ (GEOS-backed get_coordinates), NumPy 1.24+, GeoPandas 0.14+, PyArrow 14+, and SciPy 1.11+ (only for the Savitzky-Golay fallback).
• Input geometry is a shapely.geometry.LineString (or a GeoSeries of them) in a projected, metric CRS such as a UTM zone (e.g. EPSG:32633). Geographic degrees must be reprojected first — see coordinate reference systems for AVs. Curvature on degree-based coordinates is scaled by roughly 111,320 m/° and is unusable downstream.
• Upstream stage: vectorized lane centerlines produced by the centerline generation algorithms stage, persisted as GeoParquet tiles. Vertices should already be ordered along the direction of travel.
• Output: a 1-D float64 array of $\kappa$ values in m⁻¹, one per vertex, ready to attach to the lane attribute record alongside heading and cross-slope.

Step-by-Step #

1. Compute curvature for a single projected LineString #

Iterating Shapely geometries vertex-by-vertex in Python incurs per-call GEOS overhead. Bulk-extract the coordinate array with shapely.get_coordinates(), which returns a contiguous (N, 2) array, then drive everything through NumPy. np.gradient applies central differences at interior points and one-sided differences at the endpoints, using the arc-length array as the (non-uniform) spacing.

import numpy as np
from shapely import get_coordinates
from shapely.geometry import LineString

def compute_curvature(line: LineString, eps: float = 1e-9) -> np.ndarray:
    """Vectorized discrete curvature for a projected Shapely LineString.

    Args:
        line: A Shapely 2.0 LineString in a metric (projected) CRS.
        eps:  Stability guard for near-zero denominator (straight segments).

    Returns:
        1-D array of κ values (m⁻¹), same length as the vertex count.
    """
    coords = get_coordinates(line)          # (N, 2), no Python loop
    dx = np.diff(coords[:, 0])
    dy = np.diff(coords[:, 1])
    ds_seg = np.hypot(dx, dy)               # segment lengths
    s = np.concatenate([[0.0], np.cumsum(ds_seg)])  # cumulative arc length

    x_prime  = np.gradient(coords[:, 0], s)
    y_prime  = np.gradient(coords[:, 1], s)
    x_dprime = np.gradient(x_prime, s)
    y_dprime = np.gradient(y_prime, s)

    numerator   = np.abs(x_prime * y_dprime - y_prime * x_dprime)
    denominator = (x_prime**2 + y_prime**2) ** 1.5
    return numerator / (denominator + eps)

Key parameters: eps = 1e-9 keeps straight segments (zero denominator) from producing inf/nan; the cumulative arc length s — not the raw vertex index — must be passed to np.gradient so non-uniform vertex spacing is handled correctly. Expected output for an N-vertex line is an array of shape (N,). See the NumPy gradient reference for the spacing semantics and the Shapely 2.0 documentation for coordinate handling.

2. Stream large tile datasets in chunks #

AV spatial pipelines routinely process 50–200 GB of GeoParquet tiles. Loading every centerline at once OOMs a 32 GB worker. Read in row-group batches with pyarrow.parquet.ParquetFile.iter_batches() — not read_table().to_pandas().iter_batches(), which does not exist on a pandas DataFrame.

import pyarrow.parquet as pq
import geopandas as gpd
import numpy as np

def process_curvature_from_parquet(
    path: str,
    batch_size: int = 50_000,
) -> list[np.ndarray]:
    """Compute curvature for all centerlines in a GeoParquet file, in chunks."""
    pf = pq.ParquetFile(path)
    results = []

    for batch in pf.iter_batches(batch_size=batch_size):
        gdf = gpd.GeoDataFrame.from_arrow(batch)
        for geom in gdf.geometry:
            if geom is not None and not geom.is_empty:
                results.append(compute_curvature(geom))

    return results

Key parameter: batch_size caps the number of features resident in RAM per iteration; 50,000 features keeps peak footprint predictable on a 32 GB node. When a single polyline is split across batch boundaries, carry a 2–3 vertex overlap buffer between adjacent chunks so the gradient stays continuous and no spurious curvature spike appears at the seam.

3. Preprocess unstable geometry before differentiation #

Real lane centerlines carry collinear runs, sharp kinks, and GPS-derived jitter that destabilize finite differences. Decimate redundant vertices with shapely.simplify() at a 0.05–0.1 m tolerance, then suppress high-frequency noise with a Savitzky-Golay filter before calling compute_curvature.

from scipy.signal import savgol_filter
from shapely import get_coordinates
from shapely.geometry import LineString

def stabilize(line: LineString, tol: float = 0.08,
              window: int = 7, poly: int = 3) -> LineString:
    """Decimate, then smooth, a noisy centerline prior to curvature."""
    line = line.simplify(tol, preserve_topology=False)
    coords = get_coordinates(line)
    if len(coords) <= window:           # too short to smooth safely
        return line
    xs = savgol_filter(coords[:, 0], window, poly)
    ys = savgol_filter(coords[:, 1], window, poly)
    return LineString(np.column_stack([xs, ys]))

Key parameters: tol (0.05–0.1 m) removes redundant vertices while preserving curvature continuity; window must be odd and exceed poly, and short polylines are returned untouched to avoid the filter erroring on too-few samples.

Verification & Acceptance Criteria #

Validate the discrete output against closed-form benchmarks before wiring it into the schema. A circular arc of radius $R$ has constant $\kappa = 1/R$; accept when the per-vertex result matches within ±0.001 m⁻¹.

import numpy as np
from shapely.geometry import LineString

def _arc(radius=100.0, n=200):
    t = np.linspace(0.0, np.pi / 2, n)
    return LineString(np.column_stack([radius * np.cos(t), radius * np.sin(t)]))

kappa = compute_curvature(_arc(radius=100.0))
# Ignore the two endpoints (one-sided differences are less accurate there).
interior = kappa[2:-2]
assert np.allclose(interior, 1.0 / 100.0, atol=1e-3), interior.max()
print(f"max |κ - 1/R| = {np.abs(interior - 0.01).max():.2e}")   # < 1e-3

Acceptance gate for production tiles: interior $\kappa$ within ±0.001 m⁻¹ of the analytical arc, no nan/inf values (np.isfinite(kappa).all()), and curvature magnitude bounded by design limits — flag any vertex exceeding the tightest expected radius (e.g. $\kappa > 0.2$ m⁻¹, i.e. < 5 m radius) for review against the road curvature and superelevation mapping tolerances. Cross-check spikes against the lane topology before accepting; see lane-level topology modeling.

Common Errors & Fixes #

Curvature values ~10⁵ too small (everything reads as straight). The geometry is in a geographic CRS (degrees). np.gradient differentiates degrees, scaling $\kappa$ by ~1/111,320. Reproject to a metric CRS first:

gdf = gdf.to_crs(32633)   # UTM zone 33N — pick the zone covering your tile

RuntimeWarning: invalid value encountered in divide and nan in the output. A zero-length segment (duplicate consecutive vertices) makes s non-monotonic and the denominator zero. Drop duplicate points before differentiation:

coords = get_coordinates(line)
keep = np.concatenate([[True], np.hypot(np.diff(coords[:,0]),
                                        np.diff(coords[:,1])) > 1e-6])
line = LineString(coords[keep])

AttributeError: 'DataFrame' object has no attribute 'iter_batches'. Calling pf.read_table().to_pandas().iter_batches() — iter_batches lives on ParquetFile/Table, not on a pandas DataFrame. Call pf.iter_batches(batch_size=...) directly on the ParquetFile, as in step 2.

High-frequency oscillation reads as false curvature on GPS-derived lines. Raw trajectory jitter survives differentiation and is amplified by the second derivative. Run stabilize() (step 3) — decimate then Savitzky-Golay smooth — before compute_curvature, and re-validate against the ±0.001 m⁻¹ arc gate.

Related #

• Extracting Lane Boundaries from Point Cloud Data — the upstream vectorization that produces the centerline geometry this curvature pass differentiates.
• Automating Lane-Width Attribute Sync — the sibling attribute-population stage that hangs the computed $\kappa$ onto each lane record alongside width and heading.
• Coordinate Reference Systems for AVs — why curvature must be computed in a metric (projected) CRS, and how to pick the right UTM zone before differentiating.
• Lane-Level Topology Modeling — the topology layer against which spurious curvature spikes are cross-checked before acceptance.

Up one level: Road Curvature & Superelevation Mapping — the parent workflow that consumes per-vertex $\kappa$ to derive signed curvature and cross-slope for HD map serialization.