5

This question already has an answer here:

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. Polygon/LineString intersection

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.

marked as duplicate by Matte, whyzar, Dan C, nmtoken, Fezter Feb 7 '17 at 22:19

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    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
  • 2
    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
  • 1
    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
  • 1
    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
4

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

enter image description here

 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

enter image description here

  • 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

Not the answer you're looking for? Browse other questions tagged or ask your own question.