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),]

たとえば、輪郭を描くのではなく、gdal_polygonalize.py で行われた qgis で表示されるこれらすべての白い領域の境界を探しています。

qgisのgdal_polygonalize.py等高線図

4

1 に答える 1

5

Scikit Image を教えてくれた Joe Kington に感謝します。

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 に答える