假设我有两个不相交的多边形组/“岛屿”(想想两个不相邻县的人口普查区)。我的数据可能看起来像这样:
>>> p1=Polygon([(0,0),(10,0),(10,10),(0,10)])
>>> p2=Polygon([(10,10),(20,10),(20,20),(10,20)])
>>> p3=Polygon([(10,10),(10,20),(0,10)])
>>>
>>> p4=Polygon([(40,40),(50,40),(50,30),(40,30)])
>>> p5=Polygon([(40,40),(50,40),(50,50),(40,50)])
>>> p6=Polygon([(40,40),(40,50),(30,50)])
>>>
>>> df=gpd.GeoDataFrame(geometry=[p1,p2,p3,p4,p5,p6])
>>> df
geometry
0 POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))
1 POLYGON ((10 10, 20 10, 20 20, 10 20, 10 10))
2 POLYGON ((10 10, 10 20, 0 10, 10 10))
3 POLYGON ((40 40, 50 40, 50 30, 40 30, 40 40))
4 POLYGON ((40 40, 50 40, 50 50, 40 50, 40 40))
5 POLYGON ((40 40, 40 50, 30 50, 40 40))
>>>
>>> df.plot()
我希望每个岛内的多边形都采用一个 ID(可以是任意的)来表示它的组。例如,左下角的 3 个多边形可以有 IslandID = 1,右上角的 3 个多边形可以有 IslandID = 2。
我已经开发出一种方法来执行此操作,但我想知道这是否是最好/最有效的方法。我执行以下操作:
1) 创建一个 GeoDataFrame,其几何等于多多边形一元联合中的多边形。这给了我两个多边形,每个“岛”一个。
>>> SepIslands=gpd.GeoDataFrame(geometry=list(df.unary_union))
>>> SepIslands.plot()
2) 为每个组创建一个 ID。
>>> SepIslands['IslandID']=SepIslands.index+1
3) 将岛屿空间连接到原始多边形,因此每个多边形都有适当的岛屿 ID。
>>> Final=gpd.tools.sjoin(df, SepIslands, how='left').drop('index_right',1)
>>> Final
geometry IslandID
0 POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0)) 1
1 POLYGON ((10 10, 20 10, 20 20, 10 20, 10 10)) 1
2 POLYGON ((10 10, 10 20, 0 10, 10 10)) 1
3 POLYGON ((40 40, 50 40, 50 30, 40 30, 40 40)) 2
4 POLYGON ((40 40, 50 40, 50 50, 40 50, 40 40)) 2
5 POLYGON ((40 40, 40 50, 30 50, 40 40)) 2
这确实是最好/最有效的方法吗?
最佳答案
如果每组之间的差距相当大,另一个选择是 sklearn.cluster.DBSCAN聚类多边形的质心并将它们标记为聚类。
DBSCAN 代表基于密度的噪声应用程序空间聚类,它可以将紧密排列在一起的点分组。在我们的例子中,一个岛上的多边形将聚集在同一个集群中。
这也适用于两个以上的岛屿。
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
from sklearn.cluster import DBSCAN
# Note, EPS_DISTANCE = 20 is a magic number and it needs to be
# * smaller than the gap between any two islands
# * large enough to cluster polygons in one island in same cluster
EPS_DISTANCE = 20
MIN_SAMPLE_POLYGONS = 1
p1=Polygon([(0,0),(10,0),(10,10),(0,10)])
p2=Polygon([(10,10),(20,10),(20,20),(10,20)])
p3=Polygon([(10,10),(10,20),(0,10)])
p4=Polygon([(40,40),(50,40),(50,30),(40,30)])
p5=Polygon([(40,40),(50,40),(50,50),(40,50)])
p6=Polygon([(40,40),(40,50),(30,50)])
df = gpd.GeoDataFrame(geometry=[p1, p2, p3, p4, p5, p6])
# preparation for dbscan
df['x'] = df['geometry'].centroid.x
df['y'] = df['geometry'].centroid.y
coords = df.as_matrix(columns=['x', 'y'])
# dbscan
dbscan = DBSCAN(eps=EPS_DISTANCE, min_samples=MIN_SAMPLE_POLYGONS)
clusters = dbscan.fit(coords)
# add labels back to dataframe
labels = pd.Series(clusters.labels_).rename('IslandID')
df = pd.concat([df, labels], axis=1)
> df
geometry ... IslandID
0 POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0)) ... 0
1 POLYGON ((10 10, 20 10, 20 20, 10 20, 10 10)) ... 0
2 POLYGON ((10 10, 10 20, 0 10, 10 10)) ... 0
3 POLYGON ((40 40, 50 40, 50 30, 40 30, 40 40)) ... 1
4 POLYGON ((40 40, 50 40, 50 50, 40 50, 40 40)) ... 1
5 POLYGON ((40 40, 40 50, 30 50, 40 40)) ... 1
[6 rows x 4 columns]
关于python - 在 Geopandas/Shapely 中识别独特的多边形分组,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/33440530/