Schneller Weg der Polygon-Kreuzung mit formschönen

Question

Ich habe eine große Anzahl von Polygonen (~ 100000) und versuchen, eine intelligente Möglichkeit zur Berechnung ihrer Schnittfläche mit einem regelmäßigen Gitter Zellen zu finden.

Momentan erstelle ich die Polygone und Gitterzellen mit Hilfe von formschönen (basierend auf ihren Eckkoordinaten). Dann gehe ich mit einer einfachen For-Schleife durch jedes Polygon und vergleiche es mit benachbarten Gitterzellen.

Nur ein kleines Beispiel zur Veranschaulichung der Polygone/Gitterzellen.

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: die Gitterzellen haben die Abmessungen und den Polygonen 0.25x0.25 1x1 bei max)

Eigentlich dies für ein einzelnes Polygon/Rasterzelle Combo mit etwa 0,003 Sekunden ist recht schnell. Wenn Sie diesen Code jedoch auf einer großen Anzahl von Polygonen ausführen (jeder kann Dutzende von Gitterzellen schneiden), dauert das auf meinem Rechner etwa 15+ Minuten (bis zu 30+ Minuten, abhängig von der Anzahl der sich kreuzenden Gitterzellen), was nicht akzeptabel ist. Leider habe ich keine Ahnung, wie es möglich ist, einen Code für Polygonschnittpunkte zu schreiben, um den Überlappungsbereich zu erhalten. Hast du irgendwelche Tipps? Gibt es eine Alternative zu formschön?

Answer 1

Überlegen Sie mithilfe von Rtree, welche Gitterzellen ein Polygon schneiden kann. Auf diese Weise können Sie die for-Schleife entfernen, die mit dem Array von lat/lons verwendet wird, was wahrscheinlich der langsame Teil ist.

Struktur der Code etwas wie folgt aus:

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

Es sieht aus wie die verfügbare The Shapely User Manual ist ziemlich veraltet, aber seit 2013/2014 hat Shapely strtree.py mit der Klasse STRtree hatte. Ich habe es benutzt und es scheint gut zu funktionieren.

Hier ist ein Ausschnitt aus der docstring:

STRtree ist ein R-Baum, der die Sort-Tile-Recursive Algorithmus erstellt wird. STRtree verwendet eine Sequenz von Geometrieobjekten als Initialisierungsparameter . Nach der Initialisierung kann die Abfragemethode verwendet werden, um eine räumliche Abfrage über diese Objekte zu erstellen.

>>> 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