Voie plus rapide de l'intersection de polygone avec shapely
j'ai un grand nombre de polygones (~100000) et essayez de trouver une façon intelligente de calcul de leur zone d'intersection avec une grille de cellules.
actuellement, je crée les polygones et les cellules de la grille à l'aide de shapely (basé sur leurs coordonnées de coin). Ensuite, en utilisant un simple pour-boucle je passe à travers chaque polygone et le compare aux cellules de grille voisines.
Juste un petit exemple pour illustrer les polygones/cellules de la grille.
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
(BTW: les mailles ont les dimensions 0.25x0.25 et les polygones 1x1 au max)
en fait c'est assez rapide pour un combo polygone/maille individuel avec environ 0.003 secondes. Cependant, l'exécution de ce code sur une énorme quantité de polygones (chacun pourrait croiser des douzaines de cellules de grille) prend environ 15+ minutes (jusqu'à 30+ min selon le nombre de cellules de grille qui se croisent) sur ma machine ce qui n'est pas acceptable. Malheureusement, je n'ai aucune idée de comment il est possible d'écrire un code pour l'intersection polygone obtenir la zone de chevauchement. Avez-vous des conseils? Est-il une alternative à galbées?
2 réponses
envisagez d'utiliser Rtree pour aider à identifier les mailles qu'un polygone peut croiser. De cette façon, vous pouvez supprimer la boucle for utilisée avec le tableau de lat/lons, qui est probablement la partie lente.
structurez votre code quelque chose comme ceci:
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
il semble que le Disponible Le Manuel D'Utilisation En Forme est plutôt périmé, mais depuis 2013/2014, strtree.py avec la classe STRtree. Je l'ai utilisé et il semble bien fonctionner.
Voici un extrait de la docstring:
STRtree est un arbre R qui est créé en utilisant le tri-Tile-Recursive algorithme. STRtree prend une séquence d'objets de géométrie comme initialisation paramètre. Après l'initialisation de la méthode de requête peut être utilisé pour faire un requête spatiale sur ces objets.
>>> 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