This is the first 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 uses the most basic operations with no optimizations. On the test machine, it …

basic_point_poly.py

import fiona, time
from shapely.geometry import shape
from shapely.geometry import Point

searchTime = time.time()

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

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

stnPntList = []
stnCntyList = []

for station in stationCollect:
    aPoint = shape(station['geometry'])
    aPoint.stnId = station['properties']['SYS_LOC_CO']
    stnPntList.append(aPoint)

for county in countyCollect:
    countyGeom = shape(county['geometry'])
    for stnPnt in stnPntList:
        if stnPnt.intersects(countyGeom):
            stnCntyList.append(dict([(stnPnt.stnId,county['properties']['COUNTYFIPS'])]))
        
countyCollect.close()

searchTime = time.time() - searchTime

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