How to Extract an Elevation Profile Along a Line in Python

Problem statement

You have a route — a footpath, a pipeline, a proposed cable run — and you need its elevation profile: height against distance, plus the total ascent. Here the route is the Llanberis Path up Snowdon, as mapped in OpenStreetMap, over the Copernicus DEM.

The obvious approach — read the DEM at the route's own vertices — gives an answer that depends on choices nobody wrote down:

coords = list(route.coords)
elevations = [dem[rowcol(transform, x, y)] for x, y in coords]
ascent = sum(max(0, b - a) for a, b in zip(elevations, elevations[1:]))
print(f"{ascent:.0f} m of ascent")

922 m of ascent

Sample it every 10 m with nearest-neighbour and it becomes 1,043 m. Interpolate bilinearly and it is 848 m. All three are "the total ascent of this route", and they differ by 23%.

Total ascent is not a property of a route. It is a property of a route and a sampling scheme, and a profile is only reproducible if you state both.

Quick answer

Resample the line at a fixed ground interval and interpolate the DEM bilinearly:

import numpy as np
from scipy.ndimage import map_coordinates


def fractional_rowcol(transform, xs, ys):
    """Row and column as floats, measured from cell centres — what map_coordinates expects."""
    cols, rows = ~transform * (np.asarray(xs, float), np.asarray(ys, float))
    return rows - 0.5, cols - 0.5


def profile(line, dem, transform, *, step=10.0):
    """Elevation profile at a fixed ground spacing, bilinearly interpolated."""
    n = max(2, int(np.ceil(line.length / step)) + 1)
    distances = np.linspace(0, line.length, n)
    points = [line.interpolate(d) for d in distances]

    rows, cols = fractional_rowcol(transform, [p.x for p in points], [p.y for p in points])
    z = map_coordinates(dem, [rows, cols], order=1, mode="nearest")

    return distances, z


distances, z = profile(route, dem, transform, step=10.0)
ascent = np.diff(z).clip(min=0).sum()
descent = -np.diff(z).clip(max=0).sum()
print(f"{route.length / 1000:.2f} km · {z.min():.0f}-{z.max():.0f} m · "
      f"ascent {ascent:.0f} m · descent {descent:.0f} m")

6.04 km · 242-1070 m · ascent 848 m · descent 20 m

Three parameters decide the answer:

Fixed ground spacing, bilinear sampling, gradient

1. Resample the line at a fixed ground interval

A route's own vertices are wherever the digitiser put them — dense on bends, sparse on straights. Sampling at those vertices gives a profile whose resolution varies along its length.

n = max(2, int(np.ceil(line.length / step)) + 1)
distances = np.linspace(0, line.length, n)
points = [line.interpolate(d) for d in distances]

interpolate(d) returns the point at distance d along the line, so the spacing is constant regardless of vertex density. The line must be in a projected CRS for that distance to be metres.

2. Choose the spacing from the DEM, not from taste

for step in (10, 25, 50, 100, 250):
    _, z = profile(route, dem, transform, step=step)
    print(f"every {step:4} m ({len(z):4} pts): ascent {np.diff(z).clip(min=0).sum():6.0f} m")

every   10 m ( 605 pts): ascent    848 m
every   25 m ( 243 pts): ascent    841 m
every   50 m ( 122 pts): ascent    841 m
every  100 m (  62 pts): ascent    837 m
every  250 m (  26 pts): ascent    837 m

On this route spacing barely matters: 848 m at 10 m, 837 m at 250 m. The Llanberis Path climbs almost without a break, so there are few dips for coarse sampling to skip. On a route that rises and falls, coarse sampling cuts the top off every bump and loses real ascent; test your own route rather than trusting either behaviour.

Notice where it stops mattering: between 10 m and 25 m the answer moves by 7 m out of about 850. That is the point at which you are sampling more finely than the DEM's 30 m cells, so extra samples only interpolate between the same values. A spacing of about a third of the cell size captures everything the DEM contains.

3. Use bilinear interpolation

distances = np.linspace(0, route.length, int(np.ceil(route.length / 10)) + 1)
points = [route.interpolate(d) for d in distances]
for label, (rows, cols) in [
        ("rowcol", [np.asarray(v, float) for v in rowcol(transform, [p.x for p in points], [p.y for p in points])]),
        ("fractional", fractional_rowcol(transform, [p.x for p in points], [p.y for p in points]))]:
    for order, name in [(0, "nearest"), (1, "bilinear")]:
        z = map_coordinates(dem, [rows, cols], order=order, mode="nearest")
        print(f"{label:10} {name:9}: min {z.min():.1f} max {z.max():.1f} "
              f"ascent {np.diff(z).clip(min=0).sum():.0f} m")

rowcol     nearest  : min 240.7 max 1073.6 ascent 1043 m
rowcol     bilinear : min 240.7 max 1073.6 ascent 1043 m
fractional nearest  : min 240.7 max 1073.6 ascent 1043 m
fractional bilinear : min 241.9 max 1070.0 ascent 848 m

The first two lines are the trap. rasterio.transform.rowcol returns whole row and column numbers — the cell each point falls in — so every sample sits exactly on a cell, and order=1 has nothing to interpolate: "bilinear" silently equals nearest. map_coordinates needs fractional indices measured from cell centres, which is what fractional_rowcol provides.

With real fractional positions, nearest-neighbour reports more ascent. The reason is a staircase artefact: as the sample path crosses a cell boundary the value jumps discontinuously, and every jump upward is counted as ascent. Bilinear interpolation crosses smoothly, so only genuine rises are counted.

The trade is that bilinear slightly reduces the maximum — 1,070.0 m against 1,073.6 m — because it averages the summit cell with its lower neighbours. For a peak height, sample the cell value; for a profile, interpolate.

4. Decide whether to smooth, and say so

from scipy.ndimage import uniform_filter1d

for window in (1, 3, 5, 11):
    smoothed = uniform_filter1d(z, window) if window > 1 else z
    print(f"window {window:2}: ascent {np.diff(smoothed).clip(min=0).sum():6.0f} m")

window  1: ascent    848 m
window  3: ascent    841 m
window  5: ascent    838 m
window 11: ascent    829 m

Smoothing the profile removes DEM noise, which otherwise appears as a large number of tiny rises that accumulate into real ascent. Whether that noise should count is a judgement: for a cycling route it should not, for a drainage analysis it might.

Whatever you choose, put it in the output. "848 m of ascent" means nothing; "848 m sampled every 10 m, bilinear, unsmoothed" is reproducible.

5. Compute gradient along the profile

gradient = np.gradient(z, distances)               # rise per metre
grade_pct = gradient * 100

print(f"steepest climb  {grade_pct.max():5.1f}%")
print(f"steepest descent {grade_pct.min():5.1f}%")
print(f"over 10%: {(np.abs(grade_pct) > 10).mean():.1%} of the route")

steepest climb   39.7%
steepest descent -19.2%
over 10%: 62.8% of the route

Gradient is far more sensitive to sample spacing than total ascent is, because it is a derivative. A 10 m spacing on a 30 m DEM produces gradients between interpolated points rather than between real measurements — smooth before differentiating, or sample at the cell size.

A profile of the Llanberis Path, and its sensitivity

The data used on this page

The Llanberis Path from OpenStreetMap — route relation 4004181, fetched from the OSM API and merged into one line — and the Copernicus DEM (GLO-30) around it, read from its public bucket and reprojected to British National Grid at 30 m, so that distances along the route and the cell size are both in metres.

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.merge import merge
from rasterio.transform import rowcol
from rasterio.warp import calculate_default_transform, reproject, Resampling
from shapely.geometry import LineString
from shapely.ops import linemerge

osm = requests.get("https://api.openstreetmap.org/api/0.6/relation/4004181/full.json", timeout=60,
                   headers={"User-Agent": "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"}).json()
nodes = {e["id"]: (e["lon"], e["lat"]) for e in osm["elements"] if e["type"] == "node"}
ways = [LineString([nodes[n] for n in e["nodes"]]) for e in osm["elements"] if e["type"] == "way"]
gpd.GeoDataFrame({"name": ["Llanberis Path"]}, geometry=[linemerge(ways)], crs=4326).to_file("route.gpkg")

BBOX = (-4.16, 53.05, -4.05, 53.12)
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
       "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
mosaic, src_transform = merge([COP.format("W005"), COP.format("W004")], bounds=BBOX)
transform, width, height = calculate_default_transform("EPSG:4326", "EPSG:27700", mosaic.shape[2],
                                                       mosaic.shape[1], *BBOX, resolution=30)
dem = np.full((height, width), np.nan, dtype="float32")
reproject(mosaic[0], dem, src_transform=src_transform, src_crs="EPSG:4326", dst_transform=transform,
          dst_crs="EPSG:27700", dst_nodata=np.nan, resampling=Resampling.bilinear)
with rasterio.open("dem_bng.tif", "w", driver="GTiff", width=width, height=height, count=1,
                   dtype="float32", crs="EPSG:27700", transform=transform, nodata=np.nan) as dst:
    dst.write(dem, 1)
dem = dem.astype("float64")
route = gpd.read_file("route.gpkg").to_crs("EPSG:27700").geometry.iloc[0]
print(f"route {route.length / 1000:.2f} km, {len(route.coords)} vertices; "
      f"DEM {width} × {height} at 30 m")

route 6.04 km, 276 vertices; DEM 254 × 267 at 30 m

Example 1 — a complete, reproducible profile

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import rowcol
from scipy.ndimage import map_coordinates, uniform_filter1d


def elevation_profile(line, dem_path, *, step=None, order=1, smooth=1,
                      nodata_below=-100):
    """Elevation profile with every sampling decision recorded in the result."""
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype("float64")
        if src.nodata is not None:
            dem = np.where(dem == src.nodata, np.nan, dem)
        transform, crs = src.transform, src.crs
        cell = abs(src.transform.a)

    dem = np.where(dem < nodata_below, np.nan, dem)
    if np.isnan(dem).any():
        dem = np.where(np.isnan(dem), np.nanmin(dem), dem)   # keep interpolation finite

    step = step or cell / 3
    n = max(2, int(np.ceil(line.length / step)) + 1)
    distances = np.linspace(0, line.length, n)
    points = [line.interpolate(d) for d in distances]

    rows, cols = fractional_rowcol(transform, [p.x for p in points], [p.y for p in points])
    z = map_coordinates(dem, [rows, cols], order=order, mode="nearest")
    if smooth > 1:
        z = uniform_filter1d(z, smooth)

    rise = np.diff(z)
    result = {
        "distance_m": distances,
        "elevation_m": z,
        "length_m": float(line.length),
        "min_m": float(z.min()),
        "max_m": float(z.max()),
        "ascent_m": float(rise.clip(min=0).sum()),
        "descent_m": float(-rise.clip(max=0).sum()),
        "net_m": float(z[-1] - z[0]),
        # the parameters, so the numbers can be reproduced
        "sample_step_m": float(step),
        "interpolation": {0: "nearest", 1: "bilinear", 3: "cubic"}[order],
        "smooth_window": smooth,
        "dem_cell_m": float(cell),
    }
    print(f"  {result['length_m'] / 1000:.2f} km, {n} samples every {step:.1f} m "
          f"({result['interpolation']}, smooth {smooth})")
    print(f"  {result['min_m']:.0f}-{result['max_m']:.0f} m · "
          f"ascent {result['ascent_m']:.0f} m · descent {result['descent_m']:.0f} m · "
          f"net {result['net_m']:+.0f} m")
    return result


route = gpd.read_file("route.gpkg").to_crs("EPSG:27700").geometry.iloc[0]
prof = elevation_profile(route, "dem_bng.tif", step=10)

  6.04 km, 605 samples every 10.0 m (bilinear, smooth 1)
  242-1070 m · ascent 848 m · descent 20 m · net +828 m

Every parameter that affected the numbers is in the returned dictionary. That is what makes the profile a measurement rather than an opinion.

Example 2 — showing how sensitive the answer is

import pandas as pd


def ascent_sensitivity(line, dem_path, *, steps=(10, 25, 50, 100, 250),
                       orders=(0, 1)):
    rows = []
    for order in orders:
        for step in steps:
            p = elevation_profile(line, dem_path, step=step, order=order)
            rows.append({
                "interp": p["interpolation"],
                "step_m": step,
                "samples": len(p["elevation_m"]),
                "ascent_m": round(p["ascent_m"]),
            })
    frame = pd.DataFrame(rows).pivot(index="step_m", columns="interp", values="ascent_m")
    frame["inflation"] = (frame["nearest"] / frame["bilinear"] - 1).map("{:.1%}".format)
    print(frame.to_string())
    return frame


ascent_sensitivity(route, "dem_bng.tif")

  6.04 km, 605 samples every 10.0 m (nearest, smooth 1)
  241-1074 m · ascent 1043 m · descent 210 m · net +833 m
  6.04 km, 243 samples every 25.0 m (nearest, smooth 1)
  241-1074 m · ascent 954 m · descent 121 m · net +833 m
  6.04 km, 122 samples every 50.0 m (nearest, smooth 1)
  241-1074 m · ascent 859 m · descent 26 m · net +833 m
  6.04 km, 62 samples every 100.0 m (nearest, smooth 1)
  241-1074 m · ascent 841 m · descent 8 m · net +833 m
  6.04 km, 26 samples every 250.0 m (nearest, smooth 1)
  241-1074 m · ascent 838 m · descent 5 m · net +833 m
  6.04 km, 605 samples every 10.0 m (bilinear, smooth 1)
  242-1070 m · ascent 848 m · descent 20 m · net +828 m
  6.04 km, 243 samples every 25.0 m (bilinear, smooth 1)
  242-1070 m · ascent 841 m · descent 13 m · net +828 m
  6.04 km, 122 samples every 50.0 m (bilinear, smooth 1)
  242-1070 m · ascent 841 m · descent 13 m · net +828 m
  6.04 km, 62 samples every 100.0 m (bilinear, smooth 1)
  242-1070 m · ascent 837 m · descent 9 m · net +828 m
  6.04 km, 26 samples every 250.0 m (bilinear, smooth 1)
  242-1070 m · ascent 837 m · descent 9 m · net +828 m
interp  bilinear  nearest inflation
step_m                             
10           848     1043     23.0%
25           841      954     13.4%
50           841      859      2.1%
100          837      841      0.5%
250          837      838      0.1%

Across the rows, nearest-neighbour inflates the ascent — by 23% at 10 m — and the inflation all but vanishes by 250 m, because at that spacing consecutive samples are eight cells apart and almost every step is a genuine elevation change rather than a cell-boundary jump. Down the bilinear column, spacing moves the answer by only 11 m on this steadily climbing path.

The spread across this table is 837 m to 1,043 m for one route: a 25% range, entirely from sampling choices.

Example 3 — plotting the profile with the gradient

import matplotlib.pyplot as plt


def plot_profile(prof, *, grade_threshold=10.0):
    d_km = prof["distance_m"] / 1000
    z = prof["elevation_m"]
    grade = np.gradient(z, prof["distance_m"]) * 100

    fig, (ax, gx) = plt.subplots(2, 1, figsize=(10, 6), sharex=True,
                                 gridspec_kw={"height_ratios": [3, 1]})

    ax.fill_between(d_km, z.min(), z, color="#cbd5e1", alpha=0.7)
    ax.plot(d_km, z, color="#1a3a6b", lw=1.4)
    ax.set_ylabel("elevation (m)")
    ax.set_title(f"{prof['length_m'] / 1000:.2f} km · ascent {prof['ascent_m']:.0f} m · "
                 f"sampled every {prof['sample_step_m']:.0f} m, {prof['interpolation']}")

    steep = np.abs(grade) > grade_threshold
    gx.fill_between(d_km, 0, grade, where=~steep, color="#0ea5e9", alpha=0.6)
    gx.fill_between(d_km, 0, grade, where=steep, color="#ef4444", alpha=0.8)
    gx.axhline(0, color="#64748b", lw=0.8)
    gx.set_ylabel("grade (%)")
    gx.set_xlabel("distance (km)")

    print(f"  {steep.mean():.1%} of the route steeper than {grade_threshold}%")
    fig.tight_layout()
    return fig


fig = plot_profile(prof)
fig.savefig("llanberis_path_profile.png", dpi=120)

  62.8% of the route steeper than 10.0%

The sampling parameters go in the title, not a caption nobody reads. A profile chart without them cannot be compared with any other profile chart, and comparison is usually the whole point.

The two-panel layout — elevation above, gradient below — is the convention because the eye cannot read gradient off an elevation curve whose vertical scale is exaggerated, and every profile chart exaggerates the vertical scale.

Why total ascent is not well defined

Why total ascent is not well defined

Ascent is the sum of positive elevation differences between consecutive samples. Add samples and you find more small rises; remove them and small rises average out.

At the fine end, the limit is the DEM's resolution: sampling every metre on a 30 m grid adds interpolated points between the same cell values, and the total stops changing. That is what the 10 m and 25 m rows show — 848 m and 841 m, essentially identical.

At the coarse end there is no limit. Sample a mountain route every kilometre and you get the net height change plus a couple of bumps.

The same phenomenon makes coastline length ill-defined, and the resolution here is the same: fix the measurement scale, state it, and compare only like with like. A published ascent figure without a sampling interval is not reproducible.

Why nearest-neighbour inflates ascent

The sample path crosses cell boundaries. With order=0 the returned elevation is the containing cell's value, so it is constant within a cell and jumps at the boundary — a staircase.

Every upward step is counted as ascent, and there are as many steps as cell crossings. On this 6 km route across a 30 m DEM that is a couple of hundred crossings, and the accumulated artificial rise is 195 m — the difference between 1,043 m and 848 m.

Bilinear interpolation weights the four surrounding cells by distance, so the value changes continuously and crossings contribute nothing spurious. Use order=0 only when the raster is categorical, where interpolating class codes would be meaningless.

Why gradient needs more care than elevation

Gradient is a derivative, so it amplifies noise. If elevation has an uncertainty of ±4 m and you sample every 10 m, the gradient uncertainty is roughly ±80% — larger than most gradients you would want to report.

Two defences. Sample gradient at the cell size or coarser, so consecutive samples are genuinely independent measurements. Or smooth the profile before differentiating, and state the window.

Reporting "the steepest section is 40%" from a 10 m sampling of a 30 m DEM — step 5's 39.7% — is reporting interpolation noise as much as terrain.

Why the CRS must be projected

line.length and line.interpolate(d) both work in the coordinates of the geometry. On EPSG:4326 the length is in degrees, so step=10 means ten degrees — the whole route becomes two samples.

Reproject the route to a metric CRS first. The DEM can stay geographic as long as the row/column conversion uses the DEM's own transform, which it does — the two coordinate systems only need to agree at the point where you convert route coordinates into DEM row/column, so reproject the sample points into the DEM's CRS if they differ.

Projected routes, float indices and MultiLineStrings

• The route must be in a projected CRS for length and interpolate to mean metres.

• map_coordinates takes [rows, cols] as floats. Passing integers silently gives you nearest-neighbour behaviour.

• NaN propagates through interpolation. Fill voids before sampling, or the profile has holes that break the ascent sum.

• rasterio.transform.rowcol returns whole cell indices. Feed them to map_coordinates and order=1 interpolates nothing — step 3's "bilinear" was identical to nearest until the indices were fractional and centred.

• order=0 is right for categorical rasters — land cover, classification — where averaging class codes is meaningless.

• A MultiLineString has no single interpolate. Merge it with shapely.ops.linemerge first, and check the result is a single LineString.

• Sampling below a third of the cell size adds nothing. The curve has already flattened.

• Peak height and profile maximum differ. Bilinear averaging lowers summits by a few metres; read the cell value for a peak.

• Vertical datum does not matter for ascent or gradient, only for absolute heights — see vertical datums explained.

Internal links

• Digital elevation models explained — the resolution that sets the sampling floor
• How to extract raster values at point locations — the point equivalent of this operation
• Raster resampling explained — nearest versus bilinear, generalised
• Vertical datums explained — why ascent is datum-independent
• How to calculate slope and aspect from a DEM in Python — gradient over an area rather than a line
• How to measure distance accurately in Python — the horizontal half of the profile
• Coordinate precision explained — how many digits an elevation deserves
• How to calculate the shortest path along a street network — producing the route in the first place
• How to extract the part of a route between two measures — cutting a line by measure rather than by distance

FAQ

Why does my total ascent change when I resample?

Because ascent is the sum of positive differences between samples, and more samples find more small rises. It is not a property of the route alone — always state the sampling interval.

What sample spacing should I use?

About a third of the DEM cell size. Below that the answer stops changing; above it you lose real ascent.

Should I use nearest or bilinear interpolation?

Bilinear for elevation. Nearest produces a staircase whose every upward step counts as ascent — 23% inflation at 10 m sampling on the Llanberis Path. Use nearest only for categorical rasters.

Why is the profile maximum lower than the peak height?

Bilinear interpolation averages the summit cell with its lower neighbours. For a peak height, read the cell value directly.

How do I get the gradient?

np.gradient(z, distances) * 100 for percent. Sample at the cell size or smooth first — gradient amplifies DEM noise dramatically.

My profile has gaps. Why?

NaN in the DEM propagating through interpolation. Fill voids before sampling.

Does the vertical datum affect the profile?

Not the ascent, descent or gradient — those are differences, and a constant offset cancels. It does affect the absolute heights on the axis.