A faster way to intersect polygons with beautiful

Question

I have a large number of polygons (~ 100000) and try to find a reasonable way to calculate their intersecting area with regular grid cells.

I am currently creating polygons and mesh cells using beautiful ones (based on their angular coordinates). Then, using a simple for loop, I go through each polygon and compare it with neighboring grid cells.

A small example illustrating polygon / grid cells.

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: mesh cells are 0.25x0.25 and 1x1 polygons at max.)

This is actually pretty fast for a single combined polygon / mesh with a speed of about 0.003 seconds. However, running this code on a huge number of polygons (each of which can cross dozens of grid cells) takes about 15 minutes (up to 30 + min depending on the number of intersecting grid cells) on my machine, which is unacceptable. Unfortunately, I have no idea how to write code to intersect polygons to get the area of ​​overlap. Do you have any tips? Is there an alternative to beautiful?

Answer 1

Consider using Rtree to determine which mesh cells can intersect with a polygon. This way you can remove the for loop used with the lat / lons array, which is probably the slow part.

Create your code something like this:

 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

The available Shapely user guide seems to be out of date, but since 2013/2014, Shapely has strtree.py with class STRtree. I used it and it works well.

Here is a snippet from docstring:

STRtree is an R-tree that is created using the Sort-Tile-Recursive algorithm. STRtree takes a sequence of geometry objects as an initialization parameter. After initialization, the query method can be used to create a spatial query on these objects.

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