11

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

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，另一个问题中为100,000）并且不希望全部循环的情况很有用。

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

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)

检查另一个非常简化的 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 万个点。

4

1 回答 1

7

我会使用 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 上进行数百万次查询。

于 2015-03-16T15:20:37.817 回答