How to efficiently determine which of thousands of polygons intersect with a linestring [duplicate]

Question

I have a shapefile that contains several thousand Polygons/MultiPolygons and another shapefile that has a variable number of LineStrings. In the example below the blue Polygons are those that intersect with the red LineString and the brown Polygons are those that do not.

I am using the Shapely intersects binary predicate to check whether each LineString intersects with any one of the Polygons:

import fiona
from shapely.geometry import LineString, Polygon, MultiPolygon

# Open each layer
poly_layer = fiona.open('polygon_layer.shp', 'r')
line_layer = fiona.open('line_layer.shp', 'r')

# Create each Polygon/MultiPolygon
the_polygons = []
for poly in poly_layer:
    if poly['geometry']['type'] == 'Polygon':
        the_polygons.append((poly['properties']['poly_id'],
                             Polygon(poly['geometry']['coordinates'][0])))
    elif poly['geometry']['type'] == 'MultiPolygon':
        all_polygons = []
        for poly_coords in poly['geometry']['coordinates']:
            all_polygons.append(Polygon(poly_coords[0]))

        the_polygons.append((poly['properties']['poly_id'], 
                             MultiPolygon(all_polygons)))
    else:
        print poly['geometry']['type']
        raw_input()

# Create each LineString
the_lines = []
for line in line_layer:
    the_lines.append(LineString(line['geometry']['coordinates']))

# Check for Polygons/MultiPolygons that the LineString intersects with
covered_polygons = {}
for poly_id, poly in the_polygons:
    for line in the_lines:
        if poly.intersects(line):
            covered_polygons[poly_id] = covered_polygons.get(poly_id, 0) + 1

The above code takes a long time to run, and being new to Fiona and Shapely I am not sure that this is the most efficient way to accomplish things. So is there a better way to do what I want?

Edit: The key to speeding up the computations was to use a spatial index. The linked question does not have what I believe is the fastest solution - using GeoPandas.

Answer 1

It is a much easier if you directly compute the shapely geometries (shape()):

 import fiona
 from shapely.geometry import LineString, Polygon, MultiPolygon, shape
 # Open each layer
 poly_layer = fiona.open('polygon_layer.shp')
 line_layer = fiona.open('line_layer.shp')
 # convert to lists of shapely geometries
 the_lines = [shape(line['geometry']) for line in line_layer]
 the_polygons = [shape(poly['geometry']) for poly in poly_layer]
 # intersections with references to the original polygon layer
 for i, poly in enumerate(the_polygons):
      for line in the_lines:
           if poly.intersects(line):
                print poly_layer[i]

And if you use GeoPandas (Python 2.7.x, 3.x) with a spatial index (rtree), it is quicker

import geopandas as gp
from geopandas.tools import sjoin
lines = gp.GeoDataFrame.from_file('line_layer.shp') 
poly = gp.GeoDataFrame.from_file('polygon_layer.shp')
intersections= gp.sjoin(poly, lines, how="inner", op='intersects')
intersections