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

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'])))
elif poly['geometry']['type'] == 'MultiPolygon':
all_polygons = []
for poly_coords in poly['geometry']['coordinates']:
all_polygons.append(Polygon(poly_coords))

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.

• Are you amenable to trying a spatial index such as 'rtree', else do you have access to postgis? (Which can let you rapidly filter your polygons through spatial indexing). – songololo Feb 7 '17 at 15:52
• Your problem is that you don´t use an index. In the answer here is nearly the same example just with the added use of an index. gis.stackexchange.com/a/119397/69528 – Matte Feb 7 '17 at 16:03
• I didn't know about spatial indexing and spatial databases. These answers have given me a path to go down. – derNincompoop Feb 7 '17 at 16:10
• Btw you are making things very complicated with your layer loading, just use something like `the_polygons = [shape(g['geometry']) for g in poly_layer]` to create Shapely geometries from the Fiona records. – bugmenot123 Feb 7 '17 at 16:11
• That is definitely something I didn't catch when reading the user manual, but it is much simpler to use. – derNincompoop Feb 7 '17 at 16:26

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
`````` • I think the part that takes the most time are the nested for loops where the `intersects` check happens. Unless I am missing something your code doesn't seem to address that. I agree that the first part (reading in the shapefiles) is much cleaner and simpler than what I have. – derNincompoop Feb 7 '17 at 16:59
• Look at the solution with GeoPandas – gene Feb 7 '17 at 17:36
• Wow, I implemented the rtree spatial index solution on my own and it was taking comparable time to my original solution (110 seconds versus 114), However the geopandas solution take 11 seconds. An order of magnitude faster! – derNincompoop Feb 7 '17 at 18:00