Maximizing Code Performance for Shapely

You need to use a spatial index. Without an spatial index, you must iterate through all the geometries. With a bounding spatial index, you iterate only through the geometries which have a chance to intersect the other geometries.

Popular bounding spatial indexes in Python:

  • R-tree index (Python modules Rtree or pyrtree)
  • Quadtree index (Quadtree module)

Examples:

  • rtree examples
  • rtree python polygon index
  • Best way to find the polygons crossed by a line
  • More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc with Fiona.
  • libspatialindex (and Rtree) Fast and fun spatial indexing for bounding boxe with Fiona
  • Using Rtree Spatial Indexing With OGR, shapely and Fiona with Fiona
  • Trying a Python R-tree implementation
  • Search Nearby using Rtree

In your script:

  1. Why
    • from fiona import collectionif you import fiona and use only fiona?
    • import numpy as np if you are not using Numpy but math?
    • import shapely if you specify from shapely.geometry import *?
  2. If sa is a list, then use sa = [shape(i['geometry']) for i in for i in study_area]. If not, you need a condition; otherwise sa will be the last feature of study_area.
  3. Why pol = grassland.next() if you read the shapefile with MultiPolygon([shape(pol['geometry']) for pol in grassland])?

Some suggestions:

You use the grassland and trap variables once, so you don't need them

gl = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('land_cover_type_of_interest.shp')]) 

and

for point in fiona.open('some_points.shp'):

When it comes to optimizing code, don't guess - profile https://stackoverflow.com/questions/582336/how-can-you-profile-a-python-script.

Just looking at your code, a similar statements :

avail_area = pt_buffer.intersection(sa).area

seems like it would get called even more than the statement you identified, since it's nested in yet another loop.

Also, you might look into prepared geometry with shapely - http://toblerity.org/shapely/manual.html#prepared-geometry-operations It supports limited geometric operations, so it may not be a drop in replacement for what you're looking to do.