几何集合与多边形相交

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) 可以工作，但存在性能问题。

[编辑] 这是我的真实情况。我的沃罗诺伊多边形： 部分放大： 更多变焦：

这是我的“polygon_to_intersect”（橙色轮廓）： 进行一些缩放：

正如 @Pieter 所提到的，交集使用“联合语义”，这意味着我的 Voronoi 多边形首先被联合起来，然后与我的 Polygon_to_Intersect 相交。

但这不是我想要的。我想要这个（结果文件中的每个 Voronoi 多边形；其中一些在 Polygon_to_Intersect 边界处相交）：

所以我使用了这段代码：

output_shapefile = 'intersected.shp'
n = 0
w = 0
i = 0
skipped = 0

with shapefile.Writer(output_shapefile) as writer:
    writer.field("ID", "N")
    for poly in myVoronoi:
        n += 1
        if not (n % 500):
            print(f'Processed {n}, within = {w}, intersected = {i}, skipped = {skipped}')
        # If the Voronoi polygon is fully inside polygon_to_intersect, don't waste time intersecting
        # Just write the polygon into the output file
        if within(poly, polygon_to_intersect):
            writer.shape(poly)
            writer.record(n)
            w += 1
        else:
            # Perform the intersection
            intersected = poly.intersection(polygon_to_intersect)
            if intersected:
                writer.shape(intersected)
                writer.record(n + 50000)
                i += 1
            else:
                skipped += 1

这工作正常，但花了将近 25 分钟才能完成。 myVoronoi 中有 18.600 多个多边形，polygon_to_intersect 有数十万个顶点。

Pieter 使用 Rtree 的建议在性能方面没有多大帮助。

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)

因为它执行的交集比我上面的代码要多得多。

所以我创建了一个更加简化的“polygon_to_intersect”：

只有 76 个顶点。

并修改了我的代码：

# Read the new simplified shapefile into memory
simplified_file = 'simplified.shp'
with shapefile.Reader(simplified_file, shapeType=5) as w:
    for shapeRec in w.shapeRecords():
        polygon_to_intersect_simplified = shape(shapeRec.shape)
        prepare(polygon_to_intersect_simplified)


output_shapefile = 'intersected.shp'
n = 0
w = 0
i = 0
skipped = 0

with shapefile.Writer(output_shapefile) as writer:
    writer.field("ID", "N")
    for poly in myVoronoi:
        n += 1
        if not (n % 500):
            print(f'Processed {n}, within = {w}, intersected = {i}, skipped = {skipped}')
        # If the Voronoi polygon is fully inside the simplified shape, don't waste time intersecting
        # Just write the polygon into the output file
        if within(poly, polygon_to_intersect_simplified):
            writer.shape(poly)
            writer.record(n)
            w += 1
        else:
            # Perform the intersection with the real (not simplified) polygon
            intersected = poly.intersection(polygon_to_intersect)
            if intersected:
                writer.shape(intersected)
                writer.record(n + 50000)
                i += 1
            else:
                skipped += 1

“within”相对于简化得多的多边形计算了 18.600 多次，速度更快。内部的多边形不相交，因此可以节省时间。 然后，那些不在“内部”的多边形与真实的 Polygon_to_Intersect 相交。

修改后的代码运行时间不到 4 分钟。

没有简化的多边形：

• 以内 = 16257
• 相交 = 2357

使用简化的多边形：

• 以内 = 13669 <--- very fast against the simplified polygon
• 相交 = 4945 <--- more intersections than the previous code but overall the execution time goes from 25 minutes to 4 minutes.

Answer 1

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

使用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)