使用 Shapely 查找多边形的最大内接矩形

Question

我试图在六个多边形内定位数百万个点。这是我的代码：

def find_shape(longitude,latitude):
    if longitude != 0 and latitude != 0:
        point = shapely.geometry.Point(longitude,latitude)
    else:
        return "Unknown"
    for current_shape in all_shapes:
        if current_shape['bounding_box'].contains(point):
            if current_shape['shape'].contains(point):
                return current_shape['properties']['ShapeName']
                break
    return "Unknown"

我已经阅读了其他有关使用 shapely 提高多边形内点查询性能的问题。他们建议使用 Rtree。但是，这似乎对于有许多多边形（36,000 in one question、100,000 in another）并且不希望全部循环的情况很有用。

如您所见，我已经设置了一个边界框。这是我的形状设置代码：

with fiona.open(SHAPEFILE) as f_col:
    all_shapes = []
    for shapefile_record in f_col:
        current_shape = {}
        current_shape['shape'] = shapely.geometry.asShape(shapefile_record['geometry'])
        minx, miny, maxx, maxy = current_shape['shape'].bounds
        current_shape['bounding_box'] = shapely.geometry.box(minx, miny, maxx, maxy)
        current_shape['properties'] = shapefile_record['properties']
        all_shapes.append(current_shape)

也检查另一个非常简化的形状是否有用，即由最大的内接矩形（或者可能是三角形）组成的形状？

检查匀称的文档，似乎没有这个功能。也许simplify() 的一些设置？当然，我总是想确保新的简化形状不会超出原始形状的范围，所以我不必在实际形状上调用contains()。我还认为我想让新的简化形状尽可能简单，以提高速度。

任何其他建议也值得赞赏。谢谢！

编辑：在等待回复时，我想到了创建满足我要求的简化形状的想法：

current_shape['simple_shape'] = current_shape['shape'].simplify(.5)
current_shape['simple_shape'] = current_shape['simple_shape'].intersection(current_shape['shape'])

这是我在测试每个点时的使用方法：

if current_shape['simple_shape'].contains(point):
    return current_shape['properties']['ShapeName']
elif current_shape['shape'].contains(point):
    return current_shape['properties']['ShapeName']

这并不完美，因为在完成必要的intersection() 之后，形状并没有那么简单。然而，这种方法已经使处理时间减少了 60%。在我的测试中，简单多边形用于 85% 的点查询。

EDIT 2：GIS StackExchange 上的另一个相关问题： Python Efficiency — Need suggestions about how to use OGR and Shapely in more efficient way。这涉及大约 3,000 个多边形中的 150 万个点。

Answer 1

我会使用 R-Tree。 但我会将你所有的点（而不是多边形的边界框）插入 R-Tree。

使用 r 树索引例如：http://toblerity.org/rtree/

from rtree import index
from rtree.index import Rtree

idx = index.Index();

// 插入一个点，即left == right && top == bottom，本质上会在索引中插入一个点条目

for current_point in all_points:
    idx.insert(current_point.id, (current_point.x, current_point.y, current_point.x, current_point.y))

// 现在循环遍历你的多边形

for current_shape in all_shapes:
   for n in idx.intersect( current_shape['bounding_box'] )
      if current_shape['shape'].contains(n):
         # now we know that n is inside the current shape

因此，您最终会在较大的 R-Tree 上进行六次查询，而不是在小型 R-Tree 上进行数百万次查询。