如何用Shapely实现两个多边形的差集操作？

现有代码可实现两个多边形的交集计算并生成交集多边形，需求与之相反：生成从大多边形中排除交集部分的带洞多边形。

已编写的交集代码如下：

from shapely.geometry import shape,Polygon,MultiPolygon,mapping
import geopandas as gpd

g1 = geomap
print(geomap)
g2 = geomap_world
print(geomap_world)

data=[]
for index, orig in g1.iterrows():
    for index2, ref in g2.iterrows():
        if ref['geometry'].intersects(orig['geometry']):                
            data.append({'geometry':ref['geometry'].intersection(orig['geometry'])}) 

df = gpd.GeoDataFrame(data,columns=['geometry'])

df结果示例：

162  POLYGON ((-2.16991 35.16840, -1.79299 34.52792...  

尝试将intersection替换为symmetric_difference，结果却等于想要排除的交集，请问问题出在哪里？

问题根源 #

symmetric_difference（对称差）的计算逻辑是取两个几何对象互不重叠的所有区域，也就是(A - B) ∪ (B - A)，它会同时保留两个多边形各自独有的部分，这和你要的“仅保留大多边形中去掉交集的部分”完全不符。

正确实现 #

你需要用Shapely的difference方法，它的作用是从当前几何对象中减去与另一个几何对象的重叠区域，正好匹配你的需求。

假设ref['geometry']是你要保留的“大多边形”，修改后的代码如下：

from shapely.geometry import shape,Polygon,MultiPolygon,mapping
import geopandas as gpd

g1 = geomap
g2 = geomap_world

data=[]
for index, orig in g1.iterrows():
    for index2, ref in g2.iterrows():
        if ref['geometry'].intersects(orig['geometry']):                
            # 从大多边形ref中减去与orig的交集部分
            data.append({'geometry': ref['geometry'].difference(orig['geometry'])}) 

df = gpd.GeoDataFrame(data,columns=['geometry'])

关键注意点 #

• 确认好谁是“大多边形”：如果你的目标大多边形是orig['geometry']，则把代码改成orig['geometry'].difference(ref['geometry'])
• 生成的带洞多边形会以Polygon类型存储在GeoDataFrame中，Shapely会自动维护多边形的洞信息；如果结果包含多个独立区域，会返回MultiPolygon类型，GeoPandas也能正常处理

内容的提问来源于stack exchange，提问作者bravopapa