Más

¿Comprender el uso de índices espaciales con RTree?

¿Comprender el uso de índices espaciales con RTree?


Tengo problemas para comprender el uso de índices espaciales con RTree.

Ejemplo: tengo 300 puntos de búfer y necesito conocer el área de intersección de cada búfer con un archivo de forma de polígono. El shapefile de polígono tiene> 20.000 polígonos. Se sugirió que usara índices espaciales para acelerar el proceso.

Entonces ... si creo un índice espacial para mi shapefile poligonal, ¿se "adjuntará" al archivo de alguna manera, o el índice será independiente? 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 a programar). Muestran cómo crear un índice leyendo puntos creados manualmente y luego consultando otros puntos creados manualmente, lo que devuelve los identificadores que están contenidos dentro de la ventana. Tiene sentido. Pero, no explican cómo se relacionaría eso con algún archivo original del que habría venido el índice.

Estoy pensando que debe ser algo como esto:

  1. Extraiga bboxes para cada entidad poligonal de mi shapefile de polígono y colóquelos en un índice espacial, dándoles una identificación que sea la misma que su identificación en el shapefile.
  2. Consulte ese índice para obtener los ID que se cruzan.
  3. Luego, vuelva a ejecutar mi intersección solo en las características en mi shapefile 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 shapefile de un punto que contiene solo una característica de punto y un shapefile de polígono que contiene> 20,000 características de polígono.

Estoy importando los shapefiles usando Fiona, agregando el índice espacial usando RTree y tratando de hacer la intersección usando Shapely.

Mi código de prueba se ve así:

# point shapefile que representa la ubicación de las trampas de estadísticas focales deseadas = fiona.open ('single_pt_speed_test.shp', 'r') #polygon shapefile que representa la cobertura terrestre de interés gl = MultiPolygon ([shape (pol ['geometry']) for pol in fiona.open ('class3_aa.shp', 'r')]) #search area areaKM2 = 20 #create empty space index idx = index.Index () #set initial search radio for buffer areaM2 = areaKM2 * 1000000 r = (matemáticas .sqrt (areaM2 / math.pi)) #crear índice espacial a partir de gl para i, forma en enumerate (gl): idx.insert (i, shape.bounds) #índice de consulta para identificadores que se cruzan con el búfer (eventualmente tendrá múltiples puntos) para el punto en las trampas: pt_buffer = forma (punto ['geometría']). buffer (r) intersect_ids = pt_buffer.intersection (idx)

Pero sigo obteniendo TypeError: el objeto 'Polygon' no se puede llamar


Esa es la esencia de la misma. El árbol R le permite hacer una primera pasada muy rápida y le brinda un conjunto de resultados que tendrán "falsos positivos" (los cuadros delimitadores pueden cruzarse cuando las geometrías no lo hacen precisamente). Luego, revisa el conjunto de candidatos (buscándolos del shapefile por su índice) y realiza una prueba de intersección matemáticamente precisa usando, por ejemplo, Shapely. Esta es la misma estrategia que se emplea en bases de datos espaciales como PostGIS.


Casi lo tienes, pero cometiste un pequeño error. Necesitas usar elintersecciónmétodo en el índice espacial, en lugar de pasar el índice alintersecciónmétodo en el punto almacenado en búfer. Una vez que haya encontrado una lista de características en las que los cuadros delimitadores se superponen, debe verificar si su punto en la zona de influencia realmente se cruza con las geometrías.

importar fiona de shapely.geometry importar mapeo importar rtree importar matemáticas areaM2 = areaKM2 * 1000000 r = (math.sqrt (areaM2 / math.pi)) # abrir ambas capas con fiona.open ('single_pt_speed_test.shp', 'r') como layer_pnt: con fiona.open ('class3_aa.shp', 'r') como layer_land: # crea un objeto de índice espacial vacío index = rtree.index.Index () # completa el índice espacial para fid, característica en layer_land.items (): geometry = shape (feature ['geometry']) idx.insert (fid, geometry.bounds) for feature in layer_pnt: # amortigua el punto geometry = shape (feature ['geometry']) geometry_buffered = geometry.buffer ( r) # obtener una lista de fids donde los cuadros delimitadores se cruzan con fids = [int (i) for i in index.intersection (geometry_buffered.bounds)] # acceder a las características a las que esos fids hacen referencia para fid en fids: feature_land = layer_land [fid] geometry_land = shape (feature_land ['geometry']) # verifica que las geometrías se cruzan, no solo sus bboxs if geometry.intersects (geometry_land): print ('¡Encontré una intersección!') # haz algo ng útil aquí

Si está interesado en encontrar puntos que estén a una distancia mínima de su clase terrestre, puede usar eldistanciaen su lugar (cambie la sección correspondiente de la anterior).

para entidad en layer_pnt: geometry = shape (feature ['geometry']) # expandir los límites por r en todas las direcciones limits = [a + b * r para a, b en zip (geometry.bounds, [-1, -1, 1, 1])] # obtener lista de fids donde los cuadros delimitadores se cruzan fids = [int (i) para i en index.intersection (geometry_buffered.bounds)] para fid en fids: feature_land = layer_land [fid] geometry_land = shape (feature_land ['geometry']) # compruebe que las geometrías estén dentro de r metros si geometry.distance (geometry_land) <= r: print ('¡Encontrado una coincidencia!')

Si está tardando mucho en crear su índice espacial y va a hacer esto más de unas pocas veces, debería considerar la posibilidad de 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 rtree usando un generador, como este:

def gen (colección): para fid, característica en colección.items (): geometría = forma (característica ['geometría']) rendimiento ((fid, geometry.bounds, None)) index = rtree.index.Index (gen ( layer_land))

Sí, esa es la idea. Aquí hay un extracto de este tutorial sobre el uso de un índice espacial r-tree en Python, usando bien formada, Fiona y geopandas:

Un árbol r representa objetos individuales y sus cuadros delimitadores (la "r" es para "rectángulo") como el nivel más bajo del índice espacial. Luego agrega los objetos cercanos y los representa con su cuadro delimitador agregado en el siguiente nivel superior del índice. En niveles aún más altos, el árbol r agrega cuadros delimitadores y los representa por su cuadro delimitador, de forma iterativa, hasta que todo se anida en un cuadro delimitador de nivel superior. Para buscar, el r-tree toma un cuadro de consulta y, comenzando en el nivel superior, ve qué cuadros delimitadores (si los hay) lo cruzan. Luego expande cada cuadro delimitador que se cruza y ve cuál de los cuadros delimitadores secundarios dentro de él se cruza con el cuadro de consulta. Esto procede de forma recursiva hasta que todos los cuadros que se cruzan se buscan hasta el nivel más bajo y devuelve los objetos coincidentes del nivel más bajo.