多边形与匀称相交的更快方法 [英] Faster way of polygon intersection with shapely

问题描述

我有大量多边形(~100000)，并尝试找到一种智能方法来计算它们与规则网格单元的相交区域.

I have a large number of polygons (~100000) and try to find a smart way of calculating their intersecting area with a regular grid cells.

目前，我正在使用 shapely(基于它们的角坐标)创建多边形和网格单元.然后，使用一个简单的 for 循环遍历每个多边形并将其与附近的网格单元格进行比较.

Currently, I am creating the polygons and the grid cells using shapely (based on their corner coordinates). Then, using a simple for-loop I go through each polygon and compare it to nearby grid cells.

只是一个说明多边形/网格单元格的小例子.

Just a small example to illustrate the polygons/grid cells.

from shapely.geometry import box, Polygon
# Example polygon 
xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]]
polygon_shape = Polygon(xy)
# Example grid cell
gridcell_shape = box(129.5, -27.0, 129.75, 27.25)
# The intersection
polygon_shape.intersection(gridcell_shape).area

(顺便说一句:网格单元的尺寸为 0.25x0.25，多边形最大为 1x1)

(BTW: the grid cells have the dimensions 0.25x0.25 and the polygons 1x1 at max)

实际上，这对于单个多边形/网格单元组合来说非常快，大约需要 0.003 秒.但是，在我的机器上，在大量多边形(每个多边形可以与数十个网格单元相交)上运行此代码需要大约 15 分钟以上(最多 30 分钟以上，具体取决于相交网格单元的数量)，这是不可接受的.不幸的是，我不知道如何为多边形相交编写代码以获得重叠区域.你有什么建议吗?除了匀称还有其他选择吗?

Actually this is quite fast for an individual polygon/grid cell combo with around 0.003 seconds. However, running this code on a huge amount of polygons (each one could intersect dozens of grid cells) takes around 15+ minutes (up to 30+ min depending on the number of intersecting grid cells) on my machine which is not acceptable. Unfortunately, I have no idea how it is possible to write a code for polygon intersection to get the area of overlap. Do you have any tips? Is there an alternative to shapely?

推荐答案

考虑使用 Rtree帮助识别多边形可能与哪些网格单元格相交.这样，您可以删除与经纬度数组一起使用的 for 循环，这可能是较慢的部分.

Consider using Rtree to help identify which grid cells that a polygon may intersect. This way, you can remove the for loop used with the array of lat/lons, which is probably the slow part.

像这样构建你的代码:

from shapely.ops import cascaded_union
from rtree import index
idx = index.Index()

# Populate R-tree index with bounds of grid cells
for pos, cell in enumerate(grid_cells):
    # assuming cell is a shapely object
    idx.insert(pos, cell.bounds)

# Loop through each Shapely polygon
for poly in polygons:
    # Merge cells that have overlapping bounding boxes
    merged_cells = cascaded_union([grid_cells[pos] for pos in idx.intersection(poly.bounds)])
    # Now do actual intersection
    print(poly.intersection(merged_cells).area)

这篇关于多边形与匀称相交的更快方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！