Buscando una forma rápida de encontrar el polígono al que pertenece un punto usando Shapely

Tengo un conjunto de ~ 36,000 polígonos que representan una partición (~ condados) del país. Mi script en python recibe muchos puntos: pointId, longitud, latitud.

Para cada punto, quiero enviar de vuelta pointId, polygonId. Para cada punto, hacer un bucle en todos los polígonos y usar myPoint.inin (myPolygon) es bastante ineficiente.

Supongo que la biblioteca bien proporcionada ofrece una mejor manera de preparar el polígono para que encontrar el polígono para un punto se convierta en un camino de árbol (país, región, subregión, …)

Aquí está mi código hasta ahora:

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 ( ejemplos ) como índice de árbol R para: (1) indexar los límites de los polígonos 36k (haga esto justo después de leer el archivo json), luego (2) encuentre muy rápidamente los recuadros delimitadores de cada polígono en su punto de interés . Luego, (3) para los cuadros delimitadores de intersección de Rtree, use shapely para usar, por ejemplo, point.within(p) para hacer el análisis real de punto en polígono. Deberías ver un salto masivo de rendimiento con esta técnica.

Es genial,

Aquí está el código de muestra:

 polygons_sf = shapefile.Reader("") 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 

Gracias