Shapely中Voronoi多边形与多边形相交的结果差异及性能优化

问题背景 #

我正在使用Shapely进行地理数据相交运算，并用PyShp将结果写入shapefile，但对intersection的工作原理存在疑惑。

步骤复现 #

1. 生成Voronoi多边形 #

通过以下代码创建点集并生成Voronoi多边形：

# (1)
from shapely import LineString, MultiPoint, Point, Polygon, MultiPolygon
from shapely.geometry import shape
from shapely.constructive import voronoi_polygons
from shapely import intersection
import shapefile

output_shapefile = "/tmp/simple_test_01.shp"
points = MultiPoint([Point(1,1), Point(2,2), Point(4,3), Point(2,4), Point(2,5), Point(6,5), Point(5,4), Point(7,1)])
myVoronoi = voronoi_polygons(points)
#>>> myVoronoi
#<GEOMETRYCOLLECTION (POLYGON ((-5 -5, -5 4.625, -4.5 4.5, 0 3, 4 -1, 4 -5, -...>
n = 0
with shapefile.Writer(output_shapefile, shapeType=5) as w:
    w.field("ID", "N")
    for geom in myVoronoi.geoms:
        w.shape(geom)  
        w.record(n)
        n += 1

生成的Voronoi多边形在ArcMap中可视化后，包含多个围绕绿色标注点的多边形区域。

2. 创建目标相交多边形 #

创建用于相交的目标多边形：

# (2)
polygon_to_intersect = Polygon([Point(0,0), Point(2,8), Point(4,3.7), Point(10,4.5), Point(5,-4), Point(0,0)])
output_shapefile = "/tmp/simple_test_02.shp"
with shapefile.Writer(output_shapefile, shapeType=5) as w:
    w.field("ID", "N")
    w.shape(polygon_to_intersect)  
    w.record(1)

该多边形在ArcMap中显示为红色轮廓的不规则多边形。

3. 直接执行整体相交（不符合预期） #

执行整体相交运算：

# (3)
intersected = intersection(myVoronoi, polygon_to_intersect)

>>> intersected
<POLYGON ((0.6 2.4, 0.75 3, 1.125 4.5, 2 8, 3.553 4.66, 3.739 4.261, 4 3.7, ...>

output_shapefile = "/tmp/simple_test_03.shp"
n = 0
with shapefile.Writer(output_shapefile) as w:
    w.field("ID", "N")
    w.shape(intersected)  
    w.record(n)

得到的结果是单个Polygon，而非预期的MultiPolygon，仅为目标多边形与所有Voronoi多边形合并后的相交区域，丢失了单个Voronoi多边形的独立相交结果。

4. 逐个Voronoi多边形相交（符合预期但性能差） #

逐个对每个Voronoi多边形执行相交运算：

# (4)
output_shapefile = "/tmp/simple_test_04.shp"
multipol = []
for geom in myVoronoi.geoms:
    multipol.append(intersection(geom, polygon_to_intersect))

n = 0
with shapefile.Writer(output_shapefile) as w:
    w.field("ID", "N")
    for pol in multipol:
        w.shape(pol)  
        w.record(n)
        n += 1

得到了预期结果：每个着色区域对应单个Voronoi多边形与目标多边形的相交部分。但在真实场景中，当Voronoi多边形数量达18000+时，该方法耗时近25分钟，存在性能问题。

5. 转换为MultiPolygon后相交（仍不符合预期） #

尝试将Voronoi的GeometryCollection转换为MultiPolygon后再执行相交：

# (5)
output_shapefile = "/tmp/simple_test_05.shp"

myVoronoi2 = MultiPolygon(myVoronoi)
>>> myVoronoi2
<MULTIPOLYGON (((-5 -5, -5 4.625, -4.5 4.5, 0 3, 4 -1, 4 -5, -5 -5)), ((-5 1...>

intersected2 = intersection(myVoronoi2, polygon_to_intersect)

>>> intersected2
<POLYGON ((0.6 2.4, 0.75 3, 1.125 4.5, 2 8, 3.553 4.66, 3.739 4.261, 4 3.7, ...>

n = 0
with shapefile.Writer(output_shapefile) as w:
    w.field("ID", "N")
    w.shape(intersected2)  
    w.record(n)

>>> intersected2 == intersected
True

结果仍与直接整体相交一致，还是单个Polygon，无法保留独立的相交结果。

原理理解 #

经提示，Shapely的intersection方法采用并集语义：当传入集合类几何（如GeometryCollection、MultiPolygon）时，会先将所有子几何合并为一个整体，再与目标几何执行相交运算，因此无法保留单个子几何的独立相交结果。

优化方案 #

尝试使用Rtree进行空间索引优化但效果不佳，后续通过以下策略将运行时间缩短至4分钟：

• 创建目标多边形的简化版本（降低几何复杂度）
• 对每个Voronoi多边形先进行快速判断：
  • 若Voronoi多边形完全在简化后的目标多边形内部，则直接保留该多边形，无需执行相交运算
  • 若不在内部，再执行与原目标多边形的相交运算

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