Skip to content

Instantly share code, notes, and snippets.

Created February 26, 2012 19:21
Show Gist options
  • Save jsundram/1918435 to your computer and use it in GitHub Desktop.
Save jsundram/1918435 to your computer and use it in GitHub Desktop.
Read a shapefile and find out which shape contains a given lat-lng
import json
from collections import defaultdict
from itertools import izip
from rtree import index
import shapefile
def read_shapefile(shape_dir):
sf = shapefile.Reader(shape_dir)
shapes = {}
for (s, sr) in izip(sf.shapes(), sf.shapeRecords()):
shapes[int(sr.record[3])] = s
return shapes
def make_index(county_shapes):
idx = index.Index()
for (countyId, s) in county_shapes.iteritems():
idx.insert(countyId, s.bbox)
return idx
def query_index(idx, shapes, lat, lng, e=.0000001):
for hit in idx.intersection((lng - e, lat - e, lng + e, lat + e)):
if hit_test(lng, lat, shapes[hit].points):
return hit
def hit_test(x, y, polygon):
inside = False
x1, y1 = polygon[0]
for i in range(len(polygon)) + [0]:
x2, y2 = polygon[i]
if min(y1, y2) < y <= max(y1, y2) and x <= max(x1, x2):
xints = (y-y1) * (x2-x1) / (y2-y1) + x1
if x1 == x2 or x <= xints:
inside = not inside
x1, y1 = x2, y2
return inside
def process(your_data):
shapes = read_shapefile('./us_county') #
idx = make_index(shapes)
summary = defaultdict(int)
for (lat, lng) in your_data:
countyId = query_index(idx, shapes, lat, lng)
summary[countyId] += 1
json.dump(summary, open('data.json', 'w')) # scale this later
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment