Read a shapefile and find out which shape contains a given lat-lng

read.py

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') # http://www2.census.gov/geo/tiger/TIGER2009/tl_2009_us_county.zip
    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