寻找一个快速的方法来找到一个点属于使用Shapely的多边形 [英] Looking for a fast way to find the polygon a point belongs to using Shapely

问题描述

我有一套约36,000多个多边形，代表该国的一个分区（县）。
我的Python脚本收到很多分：pointId，longitude，latitude。

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

我想这个整齐的库提供了更好的方法准备多边形，以便找到一个点的多边形变成一个树形路径（国家，地区，子区域...）

以下是我的代码：

 导入sys 
导入os 
导入json 
导入时间
导入字符串
 import uuid 
 
 py_id = str（uuid.uuid4（））
 
 sys.stderr.write（py_id +'\\\
'）
 sys .stderr.write（'point_in_polygon.py V131130a.\\\
'）
 sys.stderr.flush（）
 $ b $ from shapely.geometry import 
 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，obj.iteritems（）中的值：
 if key == id：
 results.append（value）
 elif not isinstance（value ，basetring）：
 _find_values（id，value）
除了AttributeError：
传递
 
尝试：
为obj中的项目：
 if不是isinstance（item，basestring）：
 _find_values（id，item）
 TypeError除外：
传递
 
如果不是isinstance（obj，basestring）：
_find_values（id，obj）
返回结果
 
 
＃1-从json加载多边形
r = find（'rings'，jsonfile）
 len_r = len（r）
 
＃2  - 从json加载属性
a = find（'attributes'，jsonfile）
 
 def insidePolygon（point，json）：
x = 0 
，而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 
返回无，无
 
 
， 
 line = sys.stdin.readline（）
如果不是行：
 break 
 try：
 args，tobedropped = string.split（line，\\\
 ，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如果未找到则为0 
＃vehicleId将为< 0如果出现异常
 if（cantonId == None）：
 print\t.join（[0，，]）
 else：
打印\t.join（[str（vehicleId），str（cantonId），str（cantonName）]）
 
除ValueError：
打印\ t .join（[ -  1，，]）
 sys.stderr.write（py_id +'\\\
'）
 sys.stderr.write（'ValueError in Python script\\ \\ n'）
 sys.stderr.write（line）
 sys.stderr.flush（）
除外：
 sys.stderr.write（py_id +'\\\
' ）
 sys.stderr.write（'Python in script'\\ n'）
 sys.stderr.write（str（sys.exc_info（）[0]）+'\\\
'）
 sys.stderr.flush（）
 print\t.join（[ -  2，，]）
  

解决方案

使用 Rtree （ R-tree索引）索引36k多边形的边界（在阅读jsonfile之后执行此操作），然后（2）非常快地找到每个多边形的相交边界框到您的兴趣点。然后，（3）对于来自Rtree的相交边界框，使用形状来使用，例如， point.within（p）来做实际的多边形点分析。你应该看到这种技术的巨大飞跃。

I have a set of ~36,000 polygons which represent a partition (~counties) of the country. My python script receives a lot of points: pointId, longitude, latitude.

For each point, I want to send back pointId, polygonId. For each point, looping into all the polygons and use myPoint.within(myPolygon) is quite inefficient.

I suppose the shapely library offers a better way to prepare the polygon so that finding the polygon for a point becomes a tree path (country, region, sub region, ...)

Here is my code so far:

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

解决方案

Use Rtree (examples) as R-tree index to: (1) index the bounds of the 36k polygons (do this just after reading jsonfile), then (2) very quickly find the intersecting bounding boxes of each polygon to your point of interest. Then, (3) for the intersecting bounding boxes from Rtree, use shapely to use, e.g. point.within(p) to do the actual point-in-polygon analysis. You should see a massive leap of performance with this technique.

这篇关于寻找一个快速的方法来找到一个点属于使用Shapely的多边形的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！