python - guia - qgis español

Manera más rápida de intersección de polígonos con bien formado

Question

Tengo una gran cantidad de polígonos (~ 100000) y trato de encontrar una manera inteligente de calcular su área de intersección con celdas de cuadrícula regulares.

Actualmente, estoy creando los polígonos y las celdas de la cuadrícula usando bien formadas (según sus coordenadas de esquina). Luego, utilizando un bucle for simple, recorro cada polígono y lo comparo con celdas de cuadrícula cercanas.

Solo un pequeño ejemplo para ilustrar los polígonos / celdas de la cuadrícula.

from shapely.geometry import box, Polygon
# Example polygon 
xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]]
polygon_shape = Polygon(xy)
# Example grid cell
gridcell_shape = box(129.5, -27.0, 129.75, 27.25)
# The intersection
polygon_shape.intersection(gridcell_shape).area

(Por cierto: las celdas de la cuadrícula tienen las dimensiones 0.25x0.25 y los polígonos 1x1 al máximo)

En realidad, esto es bastante rápido para un combo de celdas de cuadrícula / polígono individuales con alrededor de 0.003 segundos. Sin embargo, ejecutar este código en una gran cantidad de polígonos (cada uno podría intersectar docenas de celdas de cuadrícula) toma alrededor de 15 minutos (hasta 30 minutos dependiendo de la cantidad de celdas de cuadrícula que se intersectan) en mi máquina, lo cual no es aceptable. Desafortunadamente, no tengo idea de cómo es posible escribir un código para la intersección de polígonos para obtener el área de superposición. ¿Tiene algún consejo? ¿Hay una alternativa a la forma?

Answer 1

Considere usar Rtree para ayudar a identificar las celdas de la cuadrícula que un polígono puede cruzar. De esta manera, puede eliminar el bucle for utilizado con la matriz de lat / lons, que probablemente sea la parte lenta.

Estructure su código algo como esto:

from shapely.ops import cascaded_union
from rtree import index
idx = index.Index()

# Populate R-tree index with bounds of grid cells
for pos, cell in enumerate(grid_cells):
    # assuming cell is a shapely object
    idx.insert(pos, cell.bounds)

# Loop through each Shapely polygon
for poly in polygons:
    # Merge cells that have overlapping bounding boxes
    merged_cells = cascaded_union([grid_cells[pos] for pos in idx.intersection(poly.bounds)])
    # Now do actual intersection
    print poly.intersection(merged_cells).area

Answer 2

Parece que el manual de usuario disponible de The Shapely está bastante desactualizado, pero desde 2013/2014, Shapely ha tenido strtree.py con la clase STRtree. Lo he usado y parece funcionar bien.

Aquí hay un fragmento de la cadena de documentación:

STRtree es un árbol R que se crea utilizando el algoritmo de clasificación-reclasivo de mosaico. STRtree toma una secuencia de objetos de geometría como parámetro de inicialización. Después de la inicialización, el método de consulta se puede utilizar para realizar una consulta espacial sobre esos objetos.

>>> from shapely.geometry import Polygon
>>> polys = [ Polygon(((0, 0), (1, 0), (1, 1))), Polygon(((0, 1), (0, 0), (1, 0))), Polygon(((100, 100), (101, 100), (101, 101))) ]
>>> s = STRtree(polys)
>>> query_geom = Polygon(((-1, -1), (2, 0), (2, 2), (-1, 2)))
>>> result = s.query(query_geom)
>>> polys[0] in result
True