Szybszy sposób przecięcia poligonu z zgrabnym

Question

Mam dużą liczbę wielokątów (~ 100000) i staram się znaleźć inteligentny sposób obliczania ich przecinających się obszarów za pomocą zwykłych komórek siatki.

Obecnie tworzę wielokąty i komórki siatki za pomocą kształtów (na podstawie ich współrzędnych kątowych). Następnie za pomocą prostej pętli for przejdę przez każdy wielokąt i porówna go z pobliskimi komórkami siatki.

Po prostu mały przykład ilustrujący wielokąty/komórki siatki.

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: komórki siatki mają wymiary 0.25x0.25 i 1x1 wielokątów przy max)

Właściwie jest to dość szybko do indywidualnego wielokąt/combo sieci komórkowej z około 0,003 sekundy. Jednak uruchomienie tego kodu na ogromnej liczbie wielokątów (każdy z nich może przecinać dziesiątki komórek siatki) zajmuje około 15+ minut (do 30+ min w zależności od liczby przecinających się komórek siatki) na moim komputerze, co jest niedopuszczalne. Niestety, nie mam pojęcia, jak można napisać kod dla przecięcia wieloboków, aby uzyskać obszar pokrywania się. Czy masz jakieś wskazówki? Czy istnieje alternatywa dla zgrabnej?

Answer 1

Rozważ skorzystanie z Rtree, aby pomóc w zidentyfikowaniu, które komórki siatki mogą przecinać się z wielokąta. W ten sposób możesz usunąć pętlę for używaną z tablicą lat/lons, która jest prawdopodobnie wolną częścią.

Struktura Twój kod coś takiego:

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

Wygląda dostępnej The Shapely User Manual jest raczej nieaktualne, ale ponieważ 2013/2014, zgrabna była strtree.py z STRtree klasy. Użyłem go i wygląda na to, że działa dobrze.

Oto urywek z docstring:

STRtree jest R-drzewo, które są tworzone przy użyciu algorytm Sort-dachówka rekurencyjne. STRtree pobiera sekwencję obiektów geometrii jako parametr inicjalizacji . Po inicjalizacji metoda zapytania może zostać użyta do utworzenia kwerendy przestrzennej nad tymi obiektami.

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