4

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?

Improve this question
3
  • 3
    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.
    – Spacedman
    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.
    – user2950747
    Jan 17, 2018 at 18:02
  • 2
    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).
    – sgillies
    Jan 17, 2018 at 22:46

2 Answers 2

Reset to default
8

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...

enter image description here

Step 2. Apply the existing rasterio vectorize example to the quantized raster. See this example in github.

Improve this answer
3

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):
        contours = measure.find_contours(dem.read(1), height)
            for contour in contours: # find_contours returns "an ndarray of shape (n, 2), consisting of n (row, column) coordinates along the contour"
                # Convert each contour found at this height to LineString or GeoJSON and store elevation as required

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.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.