8

I have some code that I am using to determine which Shapely Polygon/MultiPolygons intersect with a number of Shapely LineStrings. Through the answers to this question the code has gone from this:

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 = [(poly['properties']['GEOID'], shape(poly['geometry'])) for poly in poly_layer]

# 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

where every possible intersection is checked, to this:

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

# 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 = [(poly['properties']['GEOID'], shape(poly['geometry'])) for poly in poly_layer]

# Create spatial index
spatial_index = rtree.index.Index()
for idx, poly_tuple in enumerate(the_polygons):
    _, poly = poly_tuple
    spatial_index.insert(idx, poly.bounds)

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

where the spatial index is used to reduce the number of intersection checks.

With the shapefiles I have (approximately 4000 polygons, and 4 lines), the original code performs 12936 .intersection() checks and takes about 114 seconds to run. The second piece of code that uses the spatial index performs only 1816 .intersection() checks but it also takes approximately 114 seconds to run.

The code to build the spatial index only takes 1-2 seconds to run, so the 1816 checks in the second piece of code are taking just about the same amount of time to perform as the 12936 checks in the original code (since the loading of shapefiles and converting to Shapely geometries is the same in both pieces of code).

I cannot see any reason that the spatial index would make the .intersects() check take longer, so I am at a loss for why this is happening.

I can only think that I am using the RTree spatial index incorrectly. Thoughts?

5

My answer is essentially based on another answer by @gene here:

More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc

He proposed the same solution using two differents methods, with or without a spatial index.

He (correctly) stated:

What is the difference ?

  • Without the index, you must iterate through all the geometries (polygons and points).
  • With a bounding spatial index (Spatial Index RTree), you iterate only through the geometries which have a chance to intersect with your current geometry ('filter' which can save a considerable amount of calculations and time...)
  • but a Spatial Index is not a magic wand. When a very large part of the dataset has to be retrieved, a Spatial Index cannot give any speed benefit.

These sentences are self-explanatory, but I preferred quoting @gene instead of proposing the same conclusions as mine (so, all the credits go to its brilliant work!).

For a better understanding of the Rtree spatial index, you may retrieve some useful informations following these links:

Another great introduction about the using of spatial indexes may be this article by @Nathan Woodrow.

  • I understand the spatial index will work best when it is able to pare down the geometries of interest to as few as possible. That is why I compared the number of geometries of interest when using the naive method (12936) to the number of geometries when using the spatial index (1816). The intersects() method takes longer when the spatial index is being used (see the time comparison above), which is why I am not sure if I am using the spatial index incorrectly. From reading the documentation and the posts linked I think I am, but I was hoping someone could point out if I was not. – derNincompoop Feb 8 '17 at 1:04
  • You final statement was: "I can only think that I am using the Rtree spatial index incorrectly.", so I thinked you was confused about the using of the spatial index and I underlineated its meaning in my answer (which wasn't off-topic). You are trying to perform some kind of statistical analysis, but the number of geometries and attempts involved should be not enough for a better undestanding of the problem. This behavior may depend from the number of geometries involved (a very little number for appreciating the power of spatial index) or from your machine. – mgri Feb 8 '17 at 9:56
4

Just to add to mgri answer.

It is important to understand what is a spatial index (How To Properly Implement a Bounding Box for Shapely & Fiona?). With my example in How to efficiently determine which of thousands of polygons intersect with a linestring

enter image description here

You can create a spatial index with the polygons

idx = index.Index()
for feat in poly_layer:
    geom = shape(feat['geometry'])
    id = int(feat['id'])
    idx.insert(id, geom.bounds,feat)

Limits of the spatial index (polygon bounds in green)

enter image description here

Or with the LineStrings

  idx = index.Index()
  for feat in line_layer:
      geom = shape(feat['geometry'])
      id = int(feat['id'])
      idx.insert(id, geom.bounds,feat)

Limits of the spatial index (LineString bound in red)

enter image description here

Now, you iterate only through the geometries which have a chance to intersect with your current geometry (in yellow)

enter image description here

I use here the LineStrings spatial index (the results are the same but with your example of 4000 polygons and 4 lines ....).

for feat1 in poly_layer:
    geom1 = shape(feat1['geometry'])
    for id in idx.intersection(geom1.bounds):
        feat2 = line_layer[id]
        geom2 = shape(feat2['geometry'])
        if geom2.intersects(geom1):
            print 'Polygon {} intersects line {}'.format(feat1['id'], feat2['id'])

  Polygon 2 intersects line 0
  Polygon 3 intersects line 0
  Polygon 6 intersects line 0
  Polygon 9 intersects line 0

You can also use a generator (example.py)

def build_ind():
     with fiona.open('polygon_layer.shp') as shapes:
         for s in shapes:
             geom = shape(s['geometry'])
             id = int(s['id'])
             yield (id, geom.bounds, s)

 p= index.Property()
 tree = index.Index(build_ind(), properties=p)
 # first line of line_layer
 first = shape(line_layer.next()['geometry'])
 # intersection of first with polygons
 tuple(tree.intersection(first.bounds))
 (6, 2, 3, 9)

You can examine the GeoPandas script sjoin.py to understand the use of Rtree.

There are many solutions but don't forget that

  • Spatial Index is not a magic wand...
  • When I use the naive method (where I perform an intersection test between every Polygon and LineString combination) I end up performing 12936 such tests. When I use the spatial index, I only have to perform 1816 tests. I believe this means that the spatial index does provide value in this use case. However, when I time the code, performing the 1816 tests takes AS LONG as performing the 12936 tests. Shouldn't the code with the spatial index be faster since there are over 11000 fewer tests being performed? – derNincompoop Feb 8 '17 at 19:39
  • So I looked into this and found that the ~11000 tests performed by only the naive code take less than 1 second to perform, while the 1816 tests performed by both sets of code take 112 seconds to perform. Now I understand what you mean by 'Spatial Index is not a magic wand' - even though it pared down the number of tests needed, the ones needed were the ones that contributed most to the time. – derNincompoop Feb 9 '17 at 1:36
0

Edit: To clarify this answer, I incorrectly believed that all the intersection tests took approximately the same amount of time. This is not the case. The reason I did not get the speed up I expected from using a spatial index is that the selection of intersection tests are the ones that took the longest to do in the first place.

As gene and mgri have already said, a spatial index is not a magic wand. Although the spatial index pared down the number of intersection tests that needed to be performed from 12936 to only 1816, the 1816 tests are the test that took the vast majority of time to compute in the first place.

The spatial index is being used correctly, but the assumption that each intersection test takes roughly the same amount of time is what is incorrect. The time required by the intersection test can vary greatly (0.05 seconds versus 0.000007 seconds).

  • 1
    You can't take into account how the spatial index influences the speed of the further intersection because it only belongs to the complexity of the geometries involved. For your case, if geometries "A" and "B" intersect in 0.05 seconds, they will still intersect in 0.05 seconds even if you previously used a spatial index (this is obviously a theoretical statement, because I think you know that the processing of anything within a processor is related to many other factors!). – mgri Feb 12 '17 at 10:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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