This is the third 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. The county features are extract…

rtree_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')

stnPntList = []
cntyPolyDict = {}
stnCntyList = []

idx = index.Index()

idxCounter = 0
for county in countyCollect:
    idx.add(idxCounter,shape(county['geometry']).bounds)
    aPoly = shape(county['geometry'])
    aPoly.cntyFips = county['properties']['COUNTYFIPS']
    aPoly.idx = idxCounter
    cntyPolyDict[idxCounter] = aPoly
    idxCounter += 1

for station in stationCollect:
    aPoint = shape(station['geometry'])
    aPoint.stnId = station['properties']['SYS_LOC_CO']
    
    hits = list(idx.intersection((aPoint.x,aPoint.y,aPoint.x,aPoint.y)))

    if len(hits) == 1:
        stnCntyList.append(dict([(aPoint.stnId,cntyPolyDict[hits[0]].cntyFips)]))

    elif len(hits) > 1:
        for hitIdx in hits:
            if cntyPolyDict[hitIdx].intersects(aPoint):
                stnCntyList.append(dict([(aPoint.stnId,cntyPolyDict[hitIdx].cntyFips)]))

searchTime = time.time() - searchTime

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