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.
The obvious approach gives an answer that depends on choices nobody wrote down:
coords = list(route.geometry.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")1412 m of ascent
Sample the same route every 250 m instead of every 10 m and it becomes 1,081 m. Use bilinear interpolation instead of nearest-neighbour and it becomes 1,203 m. All three are "the total ascent of this route", and they differ by 31%.
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 rasterio.transform import rowcol
from scipy.ndimage import map_coordinates
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 = rowcol(transform, [p.x for p in points], [p.y for p in points])
z = map_coordinates(dem, [np.asarray(rows, float), np.asarray(cols, float)],
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"{line.length / 1000:.2f} km · {z.min():.0f}-{z.max():.0f} m · "
f"ascent {ascent:.0f} m · descent {descent:.0f} m")7.98 km · 66-981 m · ascent 1203 m · descent 746 m
Three parameters decide the answer:
| Parameter | Effect | Sensible value |
|---|---|---|
| sample spacing | more samples find more ascent | one third to one half the cell size |
| interpolation | nearest inflates ascent | bilinear (order=1) |
| smoothing | removes DEM noise from the ascent | optional, state it |
Step-by-step solution
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 ( 799 pts): ascent 1203 m
every 25 m ( 320 pts): ascent 1201 m
every 50 m ( 161 pts): ascent 1192 m
every 100 m ( 81 pts): ascent 1173 m
every 250 m ( 33 pts): ascent 1081 m
Coarse sampling smooths the profile and loses ascent — 1,081 m against 1,203 m, a 10% difference from spacing alone.
Notice where it stops mattering: between 10 m and 25 m the answer moves by 2 m out of 1,200. That is the point at which you are sampling more finely than the DEM's 28 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
for order, name in [(0, "nearest"), (1, "bilinear")]:
z = map_coordinates(dem, [rows, cols], order=order, mode="nearest")
print(f"{name:9}: min {z.min():.1f} max {z.max():.1f} "
f"ascent {np.diff(z).clip(min=0).sum():.0f} m")nearest : min 65.5 max 985.6 ascent 1412 m
bilinear : min 65.5 max 981.3 ascent 1203 m
Nearest-neighbour reports 17% 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 — 981.3 m against 985.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")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. "1,202 m of ascent" means nothing; "1,202 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")Gradient is far more sensitive to sample spacing than total ascent is, because it is a derivative. A 10 m spacing on a 28 m DEM produces gradients between interpolated points rather than between real measurements — smooth before differentiating, or sample at the cell size.
Code examples
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 = rowcol(transform, [p.x for p in points], [p.y for p in points])
z = map_coordinates(dem, [np.asarray(rows, float), np.asarray(cols, float)],
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) 7.98 km, 799 samples every 10.0 m (bilinear, smooth 1)
66-981 m · ascent 1203 m · descent 746 m · net +457 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")interp bilinear nearest inflation
step_m
10 1203 1412 17.4%
25 1201 1322 10.1%
50 1192 1236 3.7%
100 1173 1198 2.1%
250 1081 1081 0.0%
Two effects are visible at once. Down the columns, coarser sampling loses ascent. Across the rows, nearest-neighbour inflates it — and the inflation vanishes entirely by 250 m, because at that spacing consecutive samples are nine cells apart and almost every step is a genuine elevation change rather than a cell-boundary jump.
The spread across this table is 1,081 m to 1,412 m for one route: a 31% 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
plot_profile(prof) 8.4% of the route steeper than 10%
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.
Explanation
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 28 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 — 1,202 m and 1,198 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 8 km route across a 28 m DEM that is roughly 290 crossings, and the accumulated artificial rise is 209 m — the difference between 1,412 m and 1,203 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 34%" from a 10 m sampling of a 28 m DEM is reporting interpolation noise.
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 rowcol 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.
Edge cases or notes
- The route must be in a projected CRS for
lengthandinterpolateto mean metres. map_coordinatestakes[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.
order=0is right for categorical rasters — land cover, classification — where averaging class codes is meaningless.- A MultiLineString has no single
interpolate. Merge it withshapely.ops.linemergefirst, 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
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 — 17% inflation at 10 m sampling in the example here. 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.