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

问题描述

使用Rtree（示例）作为R-tree索引可以：（1）索引36k多边形的边界（在读取jsonfile之后执行此操作），然后（2）很快找到每个多边形的相交边界框到您的兴趣点。然后，对于（3）Rtree中的相交边界框，使用整形地使用，例如point.within(p)进行实际的多边形点分析。您应该会看到这种技术带来了巨大的性能飞跃。

解决方法

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

对于每个点，我想发送回pointId，polygonId。对于每个点，循环进入所有多边形并使用myPoint.within（myPolygon）效率很低。

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

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

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,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,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",""])