0

I'm updating a Python script which currently calls gdal_contour as a subprocess to extract contours. I'd like to replace gdal_contour with skimage.find_contours.

However the output of the GDAL script differs slightly from scikit-image's output when used on a raster with angled sides. Here are a few images to illustrate:

Before extracting contours, my source raster is first masked (via rasterio.mask.mask) with a rectangular polygon, producing something like this:

enter image description here

When I use skimage.find_contours, it produces open contours (which is what I want) when the desired height disappears out of bounds, as in the bottom of the right edge below:

enter image description here

But when the mask is rotated, like this:

enter image description here

The contours produced by skimage.find_contours now follow the edge of the raster when they pass out of bounds and are closed, as in the bottom-right corner.

(Contours produced by GDAL from the same raster are open in the same corner.)

enter image description here

I guess I could add a margin, then clip the contours to the correct bounds, but is there another way to deal with this situation?

(I'm open to alternatives to scikit-image, see my previous question here.)

Improve this question
4
  • (I might be off here, but I think) clipping that raster with a rotated mask doesn´t mean the raster extend is clipped. a raster grid has to stay 'unrotated' and all cells outside of the mask and within the new extend (an 'upright' rectangle intersecting the corners of the rotated mask) are asssigned a (default, if not set manually) no-data value (maybe even '0'); thus, the contours will follow along the break of actual height values and the inserted no-data values, kind of treating it as a steep cliff if you will.
    – geozelot
    Apr 6, 2018 at 15:43
  • you can find out if that is the case by trying to select a pixel/grid cell just outside the mask next to that 'cliff': if it returns a value, it´s the automatically created extend. if so, and it happens to actually be '0', you could try to set the no-data value to NULL during clipping to get your open contours.
    – geozelot
    Apr 6, 2018 at 15:43
  • Thanks for the thoughts! Yes, the surrounding pixels are filled with nodata (-9999, but I've tried numpy.NaN). I think what's happening is that skimage.measure.find_contours doesn't take nodata into account (only the edge of the array) when deciding whether to open or close a contour. From the docs: "contours which intersect the array edge will be left open. All other contours will be closed." I think I might end up extracting the contours for the whole raster, then clipping by the rotated mask.
    – user2950747
    Apr 6, 2018 at 20:12
  • okay, I checked the implementation of find_contours and it seems it doesn't offer any options. consider sticking to gdal_contours since you can specify a null_value to take into account. I'm not sure if that means it returns open isohypses, though...
    – geozelot
    Apr 6, 2018 at 20:55

0

Reset to default

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.