使用Python 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"

我已经阅读了其他一些问题,这些问题涉及改善多边形点查询的性能.他们建议Rtrees.然而,这似乎是当有很多的多边形这是非常有用的情况下(36 000个问题,10万另一个),这是不可取的遍历所有.

正如您所见,我已经设置了一个边界框.这是我的形状设置代码:

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%的点查询.

编辑2:关于GIS StackExchange的另一个相关问题: Python效率 - 需要有关如何以更有效的方式使用OGR和Shapely的建议.这涉及约3,000个多边形中的150万个点.

Answer 1

我会使用R树。但是我会将所有点（而不是多边形的边界框）插入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上进行了数百万次查询。