How to extract a profile of value from a raster along a given line?

Question

How to extract a profile of values from a raster along a given shapefile line in Python?

I am struggling finding a method to extract a profile of values (e.g. topographic profile) from a raster (geotiff). The library Rasterio has a method to clip/extract value from a raster based on a polygon, but I cannot find an equivalent method for a line shapefile.

There is a basic method with scipy, but it does not inherently conserve geographic information like a method based on higher level toolbox like rasterio could provide.

In other words, I am looking for an equivalent in Python of what the tool Terrain Profile in QGIS offers.

Thanks

Chillston · Answer 1 · 2021-12-20T16:23:42.523

I had a similar problem and found a solution which works for me. The solution uses shapely to sample points on a line/lines and then accesses respective values from the GeoTiff, therefore the extracted profile follows the direction of the line. Here is the method that I ended up with:

def extract_along_line(xarr, line, n_samples=256):
    profile = []

    for i in range(n_samples):
        # get next point on the line
        point = line.interpolate(i / n_samples - 1., normalized=True)
        # access the nearest pixel in the xarray
        value = xarr.sel(x=point.x, y=point.y, method="nearest").data
        profile.append(value)
        
    return profile

Here is a working example with data from the copernicus-dem database and the line is the diagonal of the received tile:

import rioxarray
import shapely.geometry
import matplotlib.pyplot as plt

sample_tif = ('https://elevationeuwest.blob.core.windows.net/copernicus-dem/'
              'COP30_hh/Copernicus_DSM_COG_10_N35_00_E138_00_DEM.tif')

# Load xarray
tile = rioxarray.open_rasterio(sample_tif).squeeze()
# create a line (here its the diagonal of tile)
line = shapely.geometry.MultiLineString([[
            [tile.x[-1],tile.y[-1]],
            [tile.x[0], tile.y[0]]]])

# use the method from above to extract the profile
profile = extract_along_line(tile, line)
plt.plot(profile)
plt.show()

score 2 · Accepted Answer · answered Jun 10 '20 at 14:18

This is a bit different than extracting for a polygon, as you want to sample every pixel touched by the line, in the order they are touched (the polygon approaches don't care about pixel order).

It looks like it would be possible to adapt this approach to use rasterio instead. Given a line read from a shapefile using geopandas or fiona as a shapely object, you use the endpoints to derive a new equidistant projection that you use as dst_crs in a WarpedVRT and read pixel values from that. It looks like you would need to calculate the length of your line in terms of the number of pixels you want sampled, this is the width parameter of the WarpedVRT.

This approach may need to be adapted further if your line is not an approximately straight line between the endpoints.

If you want to just get the raw pixel values under the line, you should be able to use a mask in rasterio or rasterize directly, for each line. You may want to use the all_touched=True in the case of lines.

Thanks for the response. I ended up coding my own version using numpy, scipy and gdal. I create a line of points, which I then go sample with the function [get_pt_value()](https://github.com/ArcticSnow/dempy/blob/912a69f5fadb58e8496cbe302c3947ea4efc4295/dempy/rasterTool.py#L338) — Simon, Jun 11 '20 at 13:00

How to extract a profile of value from a raster along a given line?

2 Answers2