3

我正在尝试在具有较大无数据值 (1e6) 的 2D 网格中找到所有整体的边界多边形。我已经列出了使用 scipy 的标签工作的孔。在不深入 gdal 的多边形化的情况下，是否有一种简单的方法可以生成边界多边形？我看到有 matplotlib.pylab.contour，但这试图绘制一个我真的不想要的图。关于如何为每个标签获取边界多边形的任何建议（如果可能的话，最好用一种简化多边形的方法）？我敢肯定我可以写一些东西来跨越每个标记的洞的边界，但是有什么东西已经存在了吗？

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 完成的。

qgis中的gdal_polygonalize.py等高线图

4

1 回答 1

5

感谢 Joe Kington 将我指向 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)

放大左下角后：

在此处输入图像描述

于 2013-01-30T19:29:34.050 回答