Ищете быстрый способ найти многоугольник, к которому относится точка, используя Shapely
У меня есть набор из ~ 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", "", ""])