在2D网格中围绕着带有scipy标签的区域的轮廓 [英] Contours around scipy labeled regions in a 2D grid

问题描述

我试图在没有数据值(1e6)大的2D网格中找到所有整体的边界多边形.我已经得到了使用scipy标签工作的孔的清单.在不深入了解gdal的多边形的情况下，是否有一种简单的方法来生成边界多边形?我看到有matplotlib.pylab.contour，但这试图绘制一个图，我真的不想要.关于如何获取每个标签的边界多边形的任何建议(最好是如果可能的话，最好稍微简化一下多边形的方法)?我确定我可以写一些可以走过每个标记孔的边界的东西，但是已经存在了吗?

I'm trying to find the bounding polygons of all of the wholes in a 2D grid with a large no-data value (1e6). I've got the listing of holes working using scipy's label. Without dipping into gdal's polygonalize, is there an easy way to generate the bounding polygons? I see that there is matplotlib.pylab.contour, but this tries to draw a plot, which I really don't want. Any recommendation on how to get bounding polygons for each label (preferably with a way to simplify the polygons a little if possible)? I'm sure I can write something that will walk the bounds of each labelled 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's 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),]

例如而不是计算，我正在寻找qgis中所有这些白色区域的边界，这是使用gdal_polygonalize.py完成的.

e.g. rather than countouring, I'm looking for the bounds of all of these white regions as viewed in qgis, which was done with gdal_polygonalize.py.

推荐答案

感谢Joe Kington向我指出了Scikit Image.

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)

放大左下角后:

这篇关于在2D网格中围绕着带有scipy标签的区域的轮廓的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！