我有一组~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, 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", "", ""])
Mike T.. 6
使用Rtree(示例)作为R树索引:(1)索引36k多边形的边界(在读取jsonfile之后执行此操作),然后(2)非常快速地找到每个多边形的交叉边界框到您感兴趣的点.然后,(3)对于来自Rtree的交叉边界框,使用形状来使用,例如point.within(p)
进行实际的多边形点分析.你应该看到这种技术有很大的性能提升.
使用Rtree(示例)作为R树索引:(1)索引36k多边形的边界(在读取jsonfile之后执行此操作),然后(2)非常快速地找到每个多边形的交叉边界框到您感兴趣的点.然后,(3)对于来自Rtree的交叉边界框,使用形状来使用,例如point.within(p)
进行实际的多边形点分析.你应该看到这种技术有很大的性能提升.