Python 多边形与匀称相交的更快方法
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/14697442/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Faster way of polygon intersection with shapely
提问by HyperCube
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.
我有大量多边形（~100000），并尝试找到一种智能方法来计算它们与规则网格单元的相交区域。
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.
目前，我正在使用 shapely（基于它们的角坐标）创建多边形和网格单元。然后，使用一个简单的 for 循环遍历每个多边形并将其与附近的网格单元格进行比较。
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
(BTW: the grid cells have the dimensions 0.25x0.25 and the polygons 1x1 at max)
（顺便说一句：网格单元的尺寸为 0.25x0.25，多边形最大为 1x1）
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?
实际上，对于单个多边形/网格单元组合来说，这非常快，大约为 0.003 秒。但是，在我的机器上，在大量多边形（每个多边形可以与数十个网格单元相交）上运行此代码大约需要 15 分钟以上（根据相交网格单元的数量，最多需要 30 分钟以上），这是不可接受的。不幸的是，我不知道如何为多边形相交编写代码以获得重叠区域。你有什么建议吗？有没有比匀称的替代品？
采纳答案by Mike T
Consider using Rtreeto 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.
考虑使用Rtree来帮助识别多边形可能与哪些网格单元格相交。这样，您可以删除与经纬度数组一起使用的 for 循环，这可能是较慢的部分。
Structure your code something like this:
像这样构造你的代码：
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
回答by Phil
It looks like the available The Shapely User Manualis rather out of date, but since 2013/2014, Shapely has had strtree.pywith the class STRtree. I have used it and it seems to work well.
看起来可用的 The Shapely 用户手册已经过时了，但自 2013/2014 年以来，Shapely 已经有了 strtree.py和 STRtree 类。我用过它，它似乎工作得很好。
Here is a snippet from the docstring:
这是文档字符串中的一个片段：
STRtree is an R-tree that is created using the Sort-Tile-Recursive algorithm. STRtree takes a sequence of geometry objects as initialization parameter. After initialization the query method can be used to make a spatial query over those objects.
STRtree 是使用 Sort-Tile-Recursive 算法创建的 R 树。STRtree 将一系列几何对象作为初始化参数。初始化后，查询方法可用于对这些对象进行空间查询。
>>> from shapely.geometry import Polygon
>>> polys = [ Polygon(((0, 0), (1, 0), (1, 1))), Polygon(((0, 1), (0, 0), (1, 0))), Polygon(((100, 100), (101, 100), (101, 101))) ]
>>> s = STRtree(polys)
>>> query_geom = Polygon(((-1, -1), (2, 0), (2, 2), (-1, 2)))
>>> result = s.query(query_geom)
>>> polys[0] in result
True
