This is the fourth in a series of tests using fiona and shapely to read features from shapefiles and do intersect operations on them. The features represent station points in counties. There are 7142 stations and 87 counties with mildly overlapping bboxes. This example adds an RTree index based on the rtree module. In this example, the strategy …

rtree_objects_point_poly.py

import fiona, time
from shapely.geometry import shape
from shapely import speedups
from rtree import index

searchTime = time.time() 

stationPath = "C:/stations.shp"
countyPath = "C:/counties.shp"

stationCollect = fiona.open(stationPath,'r')
countyCollect = fiona.open(countyPath,'r')

stnCntyList = []

idx = index.Index()

for county in countyCollect:
    geom = shape(county['geometry'])
    idx.insert(int(county['id']), geom.bounds, obj={'properties': county['properties'], 'geom':geom})

for station in stationCollect:
    aPoint = shape(station['geometry'])
    hits = list(idx.intersection((aPoint.x,aPoint.y,aPoint.x,aPoint.y), objects="raw"))

    if len(hits) == 1:
        stnCntyList.append(dict([(station['properties']['SYS_LOC_CO'],hits[0]['properties']['COUNTYFIPS'])]))

    elif len(hits) > 1:
        for hitIdx in hits:
            if hitIdx["geom"].intersects(aPoint):
                stnCntyList.append(dict([(station['properties']['SYS_LOC_CO'],hitIdx['properties']['COUNTYFIPS'])]))

searchTime = time.time() - searchTime

print searchTime
print "Length " + str(len(stnCntyList))