# How can I extract contours from a raster with Python?

I currently run gdal_contour as a Python subprocess to extract contours from a raster file, but would like to achieve the same function with a combination of Rasterio, Shapely and Fiona.

However, once I've used Rasterio to read the raster, I'm not sure how to extract the contours.

Can anyone point me in the right direction?

• Where have you looked? A quick search found this: pypi.python.org/pypi/contours which looks a bit hacky but uses matplotlib's contour drawing code, extracting the coordinates. You are expected to do some basic searching before asking here. Jan 17, 2018 at 17:05
• Thanks for the link. I had searched and found a few ideas but not that one. Will look in more detail at it, but as you say, it does look a bit hacky. I'm trying to decide whether I should migrate from gdal_contour over to Rasterio, but won't do that unless there is a fairly straightforward approach. Jan 17, 2018 at 18:02
• Rasterio uses a few of GDAL's algorithms (the rasterizer and polygonizer), but not (yet) the contour generator (gdal.org/gdal__alg_8h.html#aceaf98ad40f159cbfb626988c054c085). Jan 17, 2018 at 22:46

I'd recommend using gdal_contour . The results are likely to be much better than any attempt to re-implement it :-)

Having said that, there should be a way to do this in rasterio. I've not tried this, so it may behave differently to how I'd expect it to work in QGIS.

Step 1. Quantize the raster into contour bands using rasterio. Use rasterio as a raster calculator. Here's the formula I use for this. This assumes a float raster. G is the elevation gap between contours. E is the existing elevation value). You might need to tweak this.

``````((int(E+G+1)/G)*G)-G
``````

That will 'posterize' a float DEM into multiples of the contour elevation (0, G, 2G, 3G...), should look like this... Step 2. Apply the existing rasterio vectorize example to the quantized raster. See this example in github.

Another approach is to use `measure.find_contours` from scikit-image.

For example, opening a raster with Rasterio and finding contours from 0 to the highest point:

``````import rasterio
from skimage import measure

with rasterio.open("dem.tif") as dem:
interval = 5
max = int(interval * round(np.amax(dem.read(1)) / interval))
for height in range(0, max, interval):
Note that you'll need to take account of the transform (`dem.transform`) when converting the contours into real-world coordinates and into a `LineString`.