Faster way of polygon intersection with shapely
Consider using Rtree to help identify which grid cells that a polygon may intersect. This way, you can remove the for loop used with the array of lat/lons, which is probably the slow part.
Structure your code something like this:
from shapely.ops import cascaded_unionfrom rtree import indexidx = index.Index()# Populate R-tree index with bounds of grid cellsfor pos, cell in enumerate(grid_cells): # assuming cell is a shapely object idx.insert(pos, cell.bounds)# Loop through each Shapely polygonfor 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)
Since 2013/2014 Shapely hasSTRtree. I have used it and it seems to work well.
Here is a snippet from the docstring:
STRtree is an R-tree that is created using the Sort-Tile-Recursivealgorithm. STRtree takes a sequence of geometry objects as initializationparameter. After initialization the query method can be used to make aspatial query over those objects.
>>> from shapely.geometry import Polygon>>> from shapely.strtree import STRtree>>> 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 resultTrue