多边形相交的更快方法 [英] 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.
当前,我正在使用shape(基于它们的角坐标)来创建多边形和网格单元.然后,使用简单的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屋!