Tengo problemas para comprender el uso de índices espaciales con RTree.
Ejemplo: tengo 300 puntos almacenados en búfer, y necesito saber el área de intersección de cada búfer con un archivo de forma poligonal. El archivo de forma poligonal tiene> 20,000 polígonos. Se sugirió que usara índices espaciales para acelerar el proceso.
SO ... Si creo un índice espacial para mi archivo de forma poligonal, ¿se "adjuntará" al archivo de alguna manera, o el índice se mantendrá solo? Es decir, después de crearlo, ¿puedo ejecutar mi función de intersección en el archivo de polígono y obtener resultados más rápidos? ¿La intersección "verá" que hay índices espaciales y sabrá qué hacer? ¿O necesito ejecutarlo en el índice y luego relacionar esos resultados con mi archivo de polígono original a través de FID o algo así?
La documentación de RTree no me está ayudando mucho (probablemente porque solo estoy aprendiendo programación). Muestran cómo crear un índice mediante la lectura de puntos creados manualmente y luego consultarlo en relación con otros puntos creados manualmente, lo que devuelve los identificadores contenidos en la ventana. Tiene sentido. Pero, no explican cómo eso se relacionaría con algún archivo original del que habría provenido el índice.
Estoy pensando que debe ser algo como esto:
- Extraiga bboxes para cada entidad poligonal del archivo de formas de mi polígono y colóquelos en un índice espacial, dándoles una identificación que sea la misma que su identificación en el archivo de formas.
- Consulta ese índice para obtener los identificadores que se cruzan.
- Luego, vuelva a ejecutar mi intersección solo en las características de mi archivo de forma original que se identificaron al consultar mi índice (no estoy seguro de cómo haría esta última parte).
¿Tengo la idea correcta? ¿Me estoy perdiendo algo?
En este momento estoy tratando de hacer que este código funcione en un archivo de forma de punto que contiene solo una entidad de punto, y un archivo de forma de polígono que contiene> 20,000 entidades de polígono.
Estoy importando los archivos de forma con Fiona, agregando el índice espacial con RTree y tratando de hacer la intersección con Shapely.
Mi código de prueba se ve así:
#point shapefile representing location of desired focal statistic
traps = fiona.open('single_pt_speed_test.shp', 'r')
#polygon shapefile representing land cover of interest
gl = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('class3_aa.shp', 'r')])
#search area
areaKM2 = 20
#create empty spatial index
idx = index.Index()
#set initial search radius for buffer
areaM2 = areaKM2 * 1000000
r = (math.sqrt(areaM2/math.pi))
#create spatial index from gl
for i, shape in enumerate(gl):
idx.insert(i, shape.bounds)
#query index for ids that intersect with buffer (will eventually have multiple points)
for point in traps:
pt_buffer = shape(point['geometry']).buffer(r)
intersect_ids = pt_buffer.intersection(idx)
Pero sigo obteniendo TypeError: el objeto 'Polígono' no es invocable
.qix
es el MapServer / GDAL / OGR / SpatiaLite árbol cuádruple índiceTypeError: 'Polygon' object is not callable
su ejemplo de actualización porque sobrescribe lashape
función que importó de forma con un objeto Polygon que crea con esta línea:for i, shape in enumerate(gl):
Respuestas:
Esa es la esencia de esto. El árbol R le permite realizar una primera pasada muy rápida y le proporciona un conjunto de resultados que tendrán "falsos positivos" (los cuadros delimitadores pueden cruzarse cuando las geometrías no lo hacen). Luego revisa el conjunto de candidatos (recuperándolos del archivo de forma por su índice) y realiza una prueba de intersección matemáticamente precisa utilizando, por ejemplo, Shapely. Esta es la misma estrategia que se emplea en bases de datos espaciales como PostGIS.
fuente
Casi lo tienes, pero has cometido un pequeño error. Debe usar el
intersection
método en el índice espacial, en lugar de pasar el índice alintersection
método en el punto almacenado. Una vez que haya encontrado una lista de entidades en las que se superponen los cuadros delimitadores, debe verificar si su punto de protección realmente intersecta las geometrías.Si está interesado en encontrar puntos que estén dentro de una distancia mínima a su clase de tierra, puede usar el
distance
método en su lugar (cambie la sección correspondiente de la anterior).Si está tomando mucho tiempo construir su índice espacial y va a hacer esto más de unas pocas veces, debería considerar serializar el índice en un archivo. La documentación describe cómo hacer esto: http://toblerity.org/rtree/tutorial.html#serializing-your-index-to-a-file
También puede ver la carga masiva de los cuadros delimitadores en el árbol utilizando un generador, como este:
fuente
Sí, esa es la idea. Aquí hay un extracto de este tutorial sobre el uso de un índice espacial r-tree en Python , usando shapely, Fiona y geopandas:
fuente