Extract an Elevation Profile Along a Line in PyQGIS

An elevation profile is a plot of height against distance along a route, and it is the question every pipeline, cycle path, sightline and cross-section starts with. QGIS has an interactive profile tool; what it does not have is a one-call algorithm that hands you the numbers. Doing it from Python takes about fifteen lines and gives you something the tool cannot: chainage, gradient, and a table you can put in a report.

This recipe belongs to Terrain & Interpolation Analysis in PyQGIS. It covers densifying the line to a sensible sampling interval, reading raster values efficiently, computing distance and gradient along the route, and turning the result into a layer or a CSV.

From a route on the map to a profileThe route is densified so that vertices sit at a regular interval along its length. Each vertex is used to read one value from the elevation raster. Plotting those values against cumulative distance produces the profile, in which a steep section on the map appears as a steep segment on the chart.The sampling interval is the resolution of the answerDEM extentvertices every 25 m along the routemchainage along the routesample no finer than the cell sizea 5 m interval on a 25 m DEM reads the same cell five timesand draws a staircase that looks like terrain

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12) or newer.
  • A DEM in a projected CRS, with nodata declared.
  • A line layer in the same CRS as the DEM. If they differ, transform the geometry first — see transforming point coordinates.

Densify the line

A route drawn with four vertices produces four samples. Densifying puts a vertex every n metres so the profile follows the ground.

import processing
from qgis.core import QgsProject

route = QgsProject.instance().mapLayersByName("proposed_route")[0]

densified = processing.run("native:densifygeometriesgivenaninterval", {
    "INPUT": route,
    "INTERVAL": 25.0,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

Breakdown: The interval is in layer units and should be at least the DEM's cell size — sampling a 25 m DEM every 5 m reads the same cell five times and draws a staircase that readers mistake for terrain. Densification adds vertices between existing ones without moving them, so the original geometry is preserved exactly and any sharp corner in the route survives. For a straight cross-section, native:pointsalonglines is the alternative and gives points directly rather than a denser line.

Sample the raster

Reading values one at a time through the data provider is the direct route and is fast enough for thousands of points.

from qgis.core import QgsPointXY

dem = QgsProject.instance().mapLayersByName("dem")[0]
provider = dem.dataProvider()
nodata = provider.sourceNoDataValue(1)

feature = next(densified.getFeatures())
vertices = feature.geometry().asPolyline()

profile = []
cumulative = 0.0
previous = None

for point in vertices:
    value, ok = provider.sample(QgsPointXY(point), 1)
    if previous is not None:
        cumulative += previous.distance(point)
    previous = point
    if ok and value != nodata:
        profile.append((cumulative, value))

Breakdown: provider.sample() returns a (value, ok) tuple, where ok is False outside the raster extent — checking it is what stops the profile silently gaining a zero where the route leaves the DEM. Comparing against the declared nodata catches gaps inside the extent, which ok does not. distance() between consecutive vertices accumulates chainage in layer units, so the first element of each tuple is metres along the route and the second is height, which is exactly the pair a profile needs.

asPolyline() assumes a single-part line. A multipart route raises or returns an empty list, so a script handling arbitrary input should call asMultiPolyline() and iterate parts, resetting or continuing the chainage depending on whether the parts are contiguous. Converting multipart to singlepart beforehand sidesteps the question entirely.

Sampling faster on long routes

Calling sample() once per vertex is fine for a few thousand points and noticeably slow for a few hundred thousand — each call opens a small window into the raster and reads it. Where a route is long or there are many routes, reading in blocks is an order of magnitude quicker.

from qgis.core import QgsRectangle

extent = feature.geometry().boundingBox()
extent.grow(dem.rasterUnitsPerPixelX() * 2)

width = int(extent.width() / dem.rasterUnitsPerPixelX())
height = int(extent.height() / dem.rasterUnitsPerPixelY())
block = provider.block(1, extent, width, height)

def value_at(point):
    col = int((point.x() - extent.xMinimum()) / dem.rasterUnitsPerPixelX())
    row = int((extent.yMaximum() - point.y()) / dem.rasterUnitsPerPixelY())
    if 0 <= col < width and 0 <= row < height:
        return block.value(row, col)
    return None

Breakdown: provider.block() reads one rectangular region into memory in a single pass, after which each lookup is arithmetic rather than I/O. The row calculation counts down from the top because raster rows are ordered north to south while y coordinates increase northwards — inverting that is the single most common bug in hand-written block indexing, and it produces a profile that is a mirror image of the real one. Growing the extent by a couple of cells avoids losing the endpoints to rounding.

This only pays off when the bounding box is a reasonable fraction of the route: a long diagonal line has a bounding box far larger than the corridor it occupies, and reading all of it into memory may cost more than the samples saved. For those, split the route into segments and read a block per segment, or stay with sample() and accept the I/O.

Compute gradient and the numbers people ask for

The profile is more useful once it carries derived values.

rows = []
for index, (chainage, height) in enumerate(profile):
    if index == 0:
        gradient = 0.0
    else:
        run = chainage - profile[index - 1][0]
        rise = height - profile[index - 1][1]
        gradient = (rise / run * 100.0) if run else 0.0
    rows.append({"chainage": chainage, "height": height, "gradient_pct": gradient})

climb = sum(r["height"] - p["height"] for p, r in zip(rows, rows[1:]) if r["height"] > p["height"])
descent = sum(p["height"] - r["height"] for p, r in zip(rows, rows[1:]) if r["height"] < p["height"])
print(f"total climb {climb:.1f} m, total descent {descent:.1f} m, "
      f"steepest {max(abs(r['gradient_pct']) for r in rows):.1f}%")

Breakdown: Gradient as a percentage is what engineering and cycling audiences expect; multiply by math.degrees(math.atan(...)) instead if the audience thinks in degrees. Total climb and descent are the summed positive and negative changes, and they are strongly dependent on the sampling interval — a finer interval picks up more noise and reports more climb, which is why two tools rarely agree on this number and why quoting the interval alongside it matters. Guarding against a zero run handles duplicate vertices, which densification occasionally produces at segment joins.

Total climb depends on how often you lookSampling the same route every hundred metres, every twenty five metres and every five metres produces increasingly detailed profiles. Each finer interval captures more small rises and falls, so the reported total climb grows even though the terrain has not changed.Always quote the sampling interval with the totalevery 100 mclimb 96 mevery 25 mclimb 118 mevery 5 mclimb 164 mon a 25 m DEM only the middle figure is supported by the data

Write the profile out

A CSV goes into a report; a point layer goes back on the map.

import csv

with open("/data/output/profile.csv", "w", newline="") as handle:
    writer = csv.DictWriter(handle, fieldnames=["chainage", "height", "gradient_pct"])
    writer.writeheader()
    writer.writerows(rows)

Breakdown: newline="" is required on Windows or every row is separated by a blank line, and it is harmless elsewhere. A CSV of chainage and height opens directly in a spreadsheet and charts in two clicks, which is usually all a report needs. Where the profile should return to the map — to symbolise the steepest sections, for instance — build a memory layer of the sampled points instead and style it with a graduated renderer on gradient, as covered in exporting an attribute table to CSV.

QGIS version compatibility

The examples target QGIS 3.34 LTR (Python 3.12).

QGIS versionPythonNotes
3.16 LTR3.7QgsRasterDataProvider.sample() and the densify algorithm present.
3.22 LTR3.9Elevation profile tool added to the interface; the API here is unchanged.
3.28 LTR3.9QgsProfileRequest and the elevation framework introduced for layer-based profiles.
3.34 LTR3.12Baseline for this page.
3.40+3.12Profile generators available for mesh and point cloud layers as well as rasters.

Troubleshooting

  • Every value is zero. The route leaves the raster extent, and the ok flag was ignored. Check it on every sample.
  • A section of the profile drops to −9999. Nodata was not compared against. Filter it out rather than plotting it.
  • The profile is a staircase. The sampling interval is finer than the cell size. Increase it to at least one cell.
  • asPolyline() returns an empty list. The feature is multipart. Use asMultiPolyline() or explode the layer first.
  • Chainage does not match the line's length. Vertices were dropped at nodata gaps but the chainage was accumulated only for kept points. Accumulate before filtering, as in the example above.
  • The values are wrong by a constant. The DEM and the route are in different CRSs and the coordinates are being read in the wrong place. Reproject one of them.

Conclusion

Densify to at least the DEM's cell size, sample through the provider while checking both the ok flag and the nodata value, accumulate chainage before filtering, and quote the sampling interval whenever you quote a total climb. The whole thing is a short loop, and the output is a table that answers questions an interactive tool cannot.

Frequently Asked Questions

Is there a built-in algorithm for this?native:pointsalonglines plus native:rastersampling gets you a point layer with values in two calls, which is often enough. The loop above exists because it gives chainage and gradient in one pass and lets you handle nodata explicitly.

How do I profile across several rasters? Build a virtual raster mosaic over them and sample that. The VRT presents a single seamless source, so the loop is unchanged — see batch reprojecting raster datasets for building consistent inputs.

Can I get a profile of something other than elevation? Yes, the code is indifferent to what the raster holds. Sampling a rainfall or noise surface along a route works identically; only the axis label changes.

Why does my total climb differ from a GPS device's? Devices apply smoothing and a minimum-change threshold before accumulating, precisely because raw accumulation is dominated by noise. Applying a small threshold — ignore changes under a metre — brings the numbers much closer together.