Python Shapely - 将 GEOMETRYCOLLECTION 与 POLYGON 相交

Question

我使用 Shapely 来交叉地理数据，并使用 Pyshp 将结果写入 shapefile。我不完全理解“交叉点”是如何工作的。

首先，我创建一些点并以这种方式计算 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

当我在 ArcMap 上绘制 shapefile 时，我得到了这个（我突出显示的绿点，而不是在 shapefile 中）：

我创建一个多边形，与 Voronoi 多边形相交：

# (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)
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)

simple_test_03.shp 是单个多边形。它的顶点可以在这里看到：

我期待“交集”能给我与 Polygon_to_Intersect 相交的 Voronoi 多边形。然而，我在多边形外部与 Voronoi 多边形相交的地方得到了与新顶点的 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, ...>

“相交”是一个多边形而不是一个多多边形，为什么？

但是，如果我这样做：

# (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

我得到了这个结果，这就是我最初想要的：

每个彩色区域都是与“polygon_to_intersect”相交的 Voronoi 多边形。

为什么我需要做（4）才能得到我想要的东西，而我期望（3）能以这种方式工作？ myVoronoi 是 GEOMETRYCOLLECTION 而不是 POLYGON 或 MULTIPOLYGON 吗？

我尝试在相交之前将 myVoronoi 转换为 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

但结果与(3)相同。

这是一个简化的示例。在我的真实情况中，“myVoronoi”有 18000 多个多边形，而“polygon_to_intersect”要大得多。 (3) 没有按我想要的方式工作，(4) 可以工作，但存在性能问题。

Answer 1

当使用 GeometryCollections 进行叠加时，内部边界首先（有意）合并在一起，因此这个结果是可以预料的。

使用 rtree 索引应该可以解决循环的性能问题，如下所示：

import numpy as np
import shapely
from shapely import MultiPoint, Polygon
points = MultiPoint([(1,1), (2,2), (4,3), (2,4), (2,5), (6,5), (5,4), (7,1)])
myVoronoi = shapely.voronoi_polygons(points)

polygon_to_intersect = Polygon([(0,0), (2,8), (4,3.7), (10,4.5), (5,-4), (0,0)])

voronoi_polys = np.array(myVoronoi.geoms)
tree = shapely.STRtree(voronoi_polys)
intersecting_idx = tree.query(polygon_to_intersect, predicate="intersects")
intersections = shapely.intersection(voronoi_polys.take(intersecting_idx), polygon_to_intersect)
print(intersections)