gdal - 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?

Answer

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.