Если вы ищете быстрый способ найти многоугольник, то точка принадлежит 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", "", ""]) 

2 Solutions collect form web for “Если вы ищете быстрый способ найти многоугольник, то точка принадлежит Shapely”

Используйте Rtree ( примеры ) в качестве индекса R-дерева для: (1) индексирования границ многоугольников 36k (сделайте это сразу после чтения jsonfile), затем (2) очень быстро найдите пересекающиеся ограничивающие прямоугольники каждого многоугольника в вашей точке интереса , Затем, (3) для пересекающихся ограничивающих прямоугольников из Rtree, используйте форму для использования, например point.within(p) чтобы выполнить фактический анализ точки в полигоне. Вы должны увидеть массивный скачок производительности с помощью этой техники.

Здорово,

Вот пример кода:

 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 

благодаря

Interesting Posts

Отсутствует таблица при запуске Django Unittest с Sqlite3

xlwt set style making error: более 4094 XF (стилей)

Неоднозначное завершение вкладки, не работающее в iPython в Windows

Как отображать китайские символы внутри рамки данных pandas?

Прозрачный фон окна (Python Tkinter)

Лучшая практика в python для возвращаемого значения при ошибке и успехе

Как очистить распечатанные заявления в IPython

Могу ли я использовать имя переменной «type» в качестве аргумента функции в Python?

сортировка многоуровневых структурированных и рекордных массивов очень медленная

Почему мне нужно дважды набирать ctrl-d?

ND-версия itertools.combinations в numpy

Непосредственно используйте библиотеку Intel mkl на Scipy разреженной матрице для вычисления точки AT с меньшим объемом памяти

Просто. * вызов * функции в python. Как исправить возврат?

«python» не распознается как внутренняя или внешняя команда

как я могу проверить подключение базы данных к mysql в django

Python - лучший язык программирования в мире.