Contours around scipy labeled areas in a 2D mesh
I am trying to find the bounding polygons of all objects in a 2D mesh with a large non-data value (1e6). I have a list of holes working using scipy label. Without going down to gdal to polygonize, is there an easy way to generate bounding polygons? I see that there is matplotlib.pylab.contour, but this is trying to make a plot that I really don't want. Any advice on how to get the bounding polygons for each label (preferably to simplify the polygons a bit)? I'm sure I can write something that will traverse the boundaries of each marked hole, but is there something that already exists?
from osgeo import gdal
from scipy import ndimage
dem_file = gdal.Open('dem.tif')
dem = dem.file.GetRasterBand(1).ReadAsArray()
# Get a binary image of the no-data regions. The no-data value is large
bin = dem > 9e5
# Find all the wholes. Anything with a label > 0.
labels, num_labels = ndimage.measurements.label(bin)
num_labels
1063
# The hole label and size. Skip 0 as that label has all the valid data.
holes = [(label, sum(labels==label)) for label in range(1, num_labels)]
holes[:3]
[(1, 7520492),
(2, 1),
(3, 1),]
eg. instead of counting, I'm looking for the boundaries of all those white areas as shown in qgis, which was done with gdal_polygonalize.py.
source to share
Thanks to Joe Kington for pointing me to Scikit Image.
from skimage import measure
contours = measure.find_contours(labels, 1)
contours[-1]
array([[ 2686.99905927, 1054. ],
[ 2686. , 1053.00094073],
[ 2685.00094073, 1054. ],
[ 2686. , 1054.99905927],
[ 2686.99905927, 1054. ]])
imshow(labels)
for n, contour in enumerate(contours):
plt.plot(contour[:,1], contour[:, 0], linewidth=2)
After enlarging in the lower left corner:
source to share