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

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 '18 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 '18 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. Apr 6 '18 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 '18 at 20:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.