Skip to content

Instantly share code, notes, and snippets.

@fawcett
Created March 27, 2013 03:06
Show Gist options
  • Save fawcett/5251267 to your computer and use it in GitHub Desktop.
Save fawcett/5251267 to your computer and use it in GitHub Desktop.
This is the second 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 the speedups and prepared geometries options. On the test machine, it to…
import fiona, time
from shapely.geometry import shape
from shapely.geometry import Point
from shapely.prepared import prep
from shapely import speedups
import sys
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 = prep(shape(county['geometry']))
for stnPnt in stnPntList:
if countyGeom.intersects(stnPnt):
stnCntyList.append(dict([(stnPnt.stnId,county['properties']['COUNTYFIPS'])]))
countyCollect.close()
searchTime = time.time() - searchTime
print searchTime
print "Length " + str(len(stnCntyList))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment