I often have to execute spatial joins between points and polygons, of say bombing events and the boundaries of the district they took place in. Quantum GIS uses ftools to execute these kinds of spatial joins, but failed on on a relatively modest join of 40k points and 9k boundaries. I could do the join in POSTGIS but I don’t want the overhead of a full spatial database for some quick analysis. The whole operation took only a few minutes to write up and few seconds to run in python, using the package pyshp to load shapefiles, rtree to build the spatial index, and shapely to do the final boundary check.
import shapefile import shapely #Load the shapefile of polygons and convert it to shapely polygon objects polygons_sf = shapefile.Reader("C:/PolygonShapeFile.shp") polygon_shapes = polygons_sf.shapes() polygon_points = [q.points for q in polygon_shapes ] from shapely.geometry import Polygon polygons = [Polygon(q) for q in polygon_points] #Load the shapefile of points and convert it to shapely point objects points_sf = shapefile.Reader("C:/PointShapeFile.shp") point_shapes = points_sf.shapes() from shapely.geometry import Point point_coords= [q.points for q in point_shapes ] points = [Point(q.points) for q in point_shapes ] #Build a spatial index based on the bounding boxes of the polygons from rtree import index idx = index.Index() count = -1 for q in polygon_shapes: count +=1 idx.insert(count, q.bbox) #Assign one or more matching polygons to each point matches =  for i in range(len(points)): #Iterate through each point temp= None print "Point ", i #Iterate only through the bounding boxes which contain the point for j in idx.intersection( point_coords[i]): #Verify that point is within the polygon itself not just the bounding box if points[i].within(polygons[j]): print "Match found! ",j temp=j break matches.append(temp) #Either the first match found, or None for no matches