在GeoDataFrame中查找不重叠的多边形

5

我有一个GeoDataFrame,其中包含一个shapely.polygons列。其中一些是不同的,一些则不是:

In [1]: gdf
Out[2]:
    geometry
1   POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))
2   POLYGON ((1 3, 1 4, 2 4, 2 3, 1 3))
3   POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))
4   POLYGON ((3 1, 3 2, 4 2, 4 1, 3 1))
5   POLYGON ((1 3, 1 4, 2 4, 2 3, 1 3))

我需要找到仅有不重叠的多边形:

In [1]: gdf_distinct
Out[2]:
    geometry
1   POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))
2   POLYGON ((1 3, 1 4, 2 4, 2 3, 1 3))
4   POLYGON ((3 1, 3 2, 4 2, 4 1, 3 1))

由于多边形不可哈希,因此我无法在Pandas中使用简单的方法:

In [1]: gdf_distinct = gdf['geometry'].unique()

TypeError: unhashable type: 'Polygon'

有没有简单有效的方法可以创建一个只包含不同多边形的新GeoDataFrame?
注:我找到了一种方法,但它只适用于完全重复的多边形,并且我认为不是非常高效。
In [1]: m = []
        for index, row in gdf.iterrows():]
           if row['geometry'] not in m:
              m.append(row['geometry'])
        gdf_distinct = GeoDataFrame(geometry=m)

首先,我尝试了df.polygon_colum.unique(),但由于多边形不可哈希,所以没有成功。现在我正在使用intersects()进行处理,当我迭代所有行以检查列中的任何多边形是否与指定的多边形相交时,这似乎不是最好的方法。 - Headshaker
1
请编辑我们的问题以包含您已尝试的代码。一个格式良好的[mcve]更有可能吸引答案。 - MechMK1
1
做了一些编辑 - 希望它能澄清情况。 - Headshaker
你想要独特的多边形还是不重叠的多边形?因为这两者并不一定相同。 - joris
在我的情况下,我想要非重叠的多边形(只是因为我的数据可能不干净) - Headshaker
如果您只想识别唯一的多边形,可以在几何图形的WKT表示上执行.unique()操作,即首先对几何图形进行序列化。 - alphabetasoup
2个回答

4

让我们从一个包含四个多边形的列表开始,其中三个多边形与其他多边形重叠:

from shapely.geometry import Polygon
import geopandas

polygons = [
    Polygon([[1, 1], [1, 3], [3, 3], [3, 1], [1, 1]]),
    Polygon([[1, 3], [1, 5], [3, 5], [3, 3], [1, 3]]),
    Polygon([[2, 2], [2, 3.5], [3.5, 3.5], [3.5, 2], [2, 2]]),
    Polygon([[3, 1], [3, 2], [4, 2], [4, 1], [3, 1]]),
]
gdf = geopandas.GeoDataFrame(data={'A': list('ABCD')}, geometry=polygons)
gdf.plot(column='A', alpha=0.75)

它们看起来像这样:

enter image description here

因此,我们可以循环遍历每个形状,然后遍历所有其他形状,并使用shapely API检查重叠。如果没有重叠,我们将其附加到输出列表中:

non_overlapping = []
for p in polygons:
    overlaps = []
    for g in filter(lambda g: not g.equals(p), polygons):
        overlaps.append(g.overlaps(p))

    if not any(overlaps):
        non_overlapping.append(p)

任何可以提供以下信息的:

['POLYGON ((3 1, 3 2, 4 2, 4 1, 3 1))']

这正是我所期望的。

但是这实际上是O(N^2),而我认为它不必如此。

因此,让我们尝试不重复检查相同的一对:

non_overlapping = []
for n, p in enumerate(polygons[:-1], 1):  # don't include the last element
    overlaps = []
    for g in polygons[n:]:  # loop from the next element to the end
        overlaps.append(g.overlaps(p))

    if not any(overlaps):
        non_overlapping.append(str(p))

我在我的电脑上得到了相同的结果,并且速度稍微更快。

我们可以通过在if语句中使用生成器而不是普通的for块来压缩循环:

non_overlapping = []
for n, p in enumerate(polygons[:-1], 1):
    if not any(p.overlaps(g) for g in polygons[n:]):
        non_overlapping.append(p)

同样的故事。

感谢您提供的出色解释和示例! - Headshaker
难道不应该是 for g in polygons[n+1:],以避免检查自身重叠吗? - alphabetasoup
Shapely的overlaps函数会在两个几何对象有一个或多个点相交时返回True,但不是所有点都相交。因此这不会产生错误的结果,但是n+1仍然会跳过n个不必要的检查。 - alphabetasoup
在我的例子中,使用n+1会将相同的多边形两次插入到“non_overlapping”列表中。这是因为我从1开始枚举,而不是从0开始。 - Paul H
啊,我错过了那个。我将其改为在GeoDataFrame上使用iterrows,所以在我的情况下n从0开始。然后我将行索引附加到列表中以供稍后的iloc使用(我需要对非重叠多边形应用一个过程,对重叠多边形应用另一个过程)。如果您在使用GeoDataFrame时有更多提示,例如使用itertuples而不是iterrows,我会很感兴趣。 - alphabetasoup
我发现这段代码是一个很好的起点,但是在应用到一个包含数百个多边形的geopandas geodataframe的真实世界示例时,我发现了一个错误:
  1. 我发现只检查从循环中的n索引到末尾的多边形(for g in polygons[n:]:)会导致重叠的多边形被包含在非重叠列表中,而非重叠的多边形则不会被包含。 我修复了这个问题,并在我的回答中添加了一些功能来回答这个问题。
- undefined

0
感谢@Paul H提供的出色答案,以及@alphabetasoup给予的深思熟虑的评论。
尽管我的解决方案并没有对问题进行不同的回答,但它是相关的。我的使用情况涉及仅查找重叠的多边形。为此,我进行了一小段代码修改,并发现需要包括最后一个元素,以免丢失其中一个重叠的多边形。
# Find polygons in a geopandas dataframe that overlap with another polygon 
# in the same dataframe as well as non-overlapping polygons
overlapping = []
non_overlapping = []
for n, p in enumerate(list(gdf.geometry)[:], 1):  # Included the last element
    overlaps = []
    for g in list(gdf.geometry)[n:]:
        overlaps.append(g.overlaps(p))
    if any(overlaps):
        overlapping.append(p)  
    if not any(overlaps):
        non_overlapping.append(p)  # Did not store as string

我的使用案例还需要保留原始geopandas geodataframe中的其他列。以下是我如何做到这一点:
overlapping = []
non_overlapping = []
for n, p in enumerate(list(gdf.geometry)[:], 0):  # Used Pythonic zero-based indexing
    if any(p.overlaps(g) for g in list(gdf.geometry)[n:]):
        # Store the index from the original dataframe
        overlapping.append(n)
    if not any(p.overlaps(g) for g in list(gdf.geometry)[n:]):
        non_overlapping.append(n)

# Create a new dataframes and reset their indexes
gdf_overlapping = gdf.iloc[overlapping]  
gdf_overlapping.reset_index(drop=True, inplace=True)
gdf_non_overlapping = gdf.iloc[non_overlapping]
gdf_non_overlapping.reset_index(drop=True, inplace=True)

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接