使用 Shapely 寻找一种快速查找点所属多边形的方法答案

【问题标题】：Looking for a fast way to find the polygon a point belongs to using Shapely使用 Shapely 寻找一种快速查找点所属多边形的方法
【发布时间】：2013-12-16 08:54:51
【问题描述】：

我有一组约 36,000 个多边形，它们代表国家的一个分区（〜县）。我的 python 脚本收到很多点：pointId、经度、纬度。

对于每个点，我想发回pointId，polygonId。对于每个点，循环到所有多边形并使用 myPoint.within(myPolygon) 效率很低。

我认为 shapely 库提供了一种更好的方法来准备多边形，以便为一个点找到多边形成为树路径（国家、地区、子地区……）

到目前为止，这是我的代码：

import sys
import os
import json
import time
import string
import uuid

py_id = str(uuid.uuid4())

sys.stderr.write(py_id + '\n')
sys.stderr.write('point_in_polygon.py V131130a.\n')
sys.stderr.flush()

from shapely.geometry import Point
from shapely.geometry import Polygon
import sys
import json
import string
import uuid
import time

jsonpath='.\cantons.json'
jsonfile = json.loads(open(jsonpath).read())

def find(id, obj):
    results = []

    def _find_values(id, obj):
        try:
            for key, value in obj.iteritems():
                if key == id:
                    results.append(value)
                elif not isinstance(value, basestring):
                    _find_values(id, value)
        except AttributeError:
            pass

        try:
            for item in obj:
                if not isinstance(item, basestring):
                    _find_values(id, item)
        except TypeError:
            pass

    if not isinstance(obj, basestring):
        _find_values(id, obj)
    return results


#1-Load polygons from json  
r=find('rings',jsonfile)
len_r = len(r)

#2-Load attributes from json
a=find('attributes',jsonfile)

def insidePolygon(point,json):
        x=0
    while x < len_r :
            y=0
            while y < len(r[x]) :
        p=Polygon(r[x][y])
        if(point.within(p)):
            return a[y]['OBJECTID'],a[y]['NAME_5']
        y=y+1
        x=x+1
    return None,None


while True:
    line = sys.stdin.readline()
    if not line:
        break
    try:
        args, tobedropped = string.split(line, "\n", 2)

        #input: vehicleId, longitude, latitude
        vehicleId, longitude, latitude = string.split(args, "\t")

        point = Point(float(longitude), float(latitude))
        cantonId,cantonName = insidePolygon(point,r)

        #output: vehicleId, cantonId, cantonName
        # vehicleId will be 0 if not found
        # vehicleId will be < 0 in case of an exception
        if (cantonId == None):
            print "\t".join(["0", "", ""])
        else:
            print "\t".join([str(vehicleId), str(cantonId), str(cantonName)])

    except ValueError:
        print "\t".join(["-1", "", ""])
        sys.stderr.write(py_id + '\n')
        sys.stderr.write('ValueError in Python script\n')
        sys.stderr.write(line)
        sys.stderr.flush()
    except:
        sys.stderr.write(py_id + '\n')
        sys.stderr.write('Exception in Python script\n')
        sys.stderr.write(str(sys.exc_info()[0]) + '\n')
        sys.stderr.flush()
        print "\t".join(["-2", "", ""])

【问题讨论】：

如果解决方案不必绑定到特定库，我只建议制作一个常规的 2D 网格，其中每个元素都是要搜索的候选多边形列表。 200x200 的网格可能是一个很好的起点，可以将要搜索的平均多边形数量减少到接近 1。
@Aswin，我在问题中添加了代码。

标签： python performance optimization graph shapely

【解决方案1】：

太好了，

这里是示例代码：

polygons_sf = shapefile.Reader("<shapefile>")
polygon_shapes = polygons_sf.shapes()
polygon_points = [q.points for q in polygon_shapes ]
polygons = [Polygon(q) for q in polygon_points]
idx = index.Index()
count = -1
for q in polygon_shapes:
    count +=1
    idx.insert(count, q.bbox)
[...]
for j in idx.intersection([point.x, point.y]):
    if(point.within(polygons[j])):
        geo1, geo2 = polygons_sf.record(j)[0], polygons_sf.record(j)[13]
        break

谢谢

【讨论】：

注意循环中的break 语句。如果您有重叠的形状，并且您想计算每个形状中的点数，这将给出较小的计数。只是为了意识到这一点而浪费了几个小时调试代码。

【解决方案2】：

使用Rtree (examples) 作为R-tree index 来： (1) 索引 36k 多边形的边界（在读取 jsonfile 后执行此操作），然后 (2) 非常快速地找到每个多边形的相交边界框多边形到您的兴趣点。然后，(3) 对于 Rtree 中的相交边界框，使用 shapely 来使用，例如point.within(p) 进行实际的多边形点分析。您应该会看到这种技术的性能大幅提升。

【讨论】：