Skip to content

Instantly share code, notes, and snippets.

@adoc
Created January 9, 2015 01:56
Show Gist options
  • Save adoc/ff558ea45aee0e9e5178 to your computer and use it in GitHub Desktop.
Save adoc/ff558ea45aee0e9e5178 to your computer and use it in GitHub Desktop.
Census data to zipcode center.
import itertools as IT
def area_of_polygon(x, y):
"""Calculates the signed area of an arbitrary polygon given its verticies
http://stackoverflow.com/a/4682656/190597 (Joe Kington)
http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm#2D%20Polygons
"""
area = 0.0
for i in range(-1, len(x) - 1):
area += x[i] * (y[i + 1] - y[i - 1])
return area / 2.0
def centroid_of_polygon(points):
"""
http://stackoverflow.com/a/14115494/190597 (mgamba)
"""
area = area_of_polygon(*zip(*points))
result_x = 0
result_y = 0
N = len(points)
points = IT.cycle(points)
x1, y1 = next(points)
for i in range(N):
x0, y0 = x1, y1
x1, y1 = next(points)
cross = (x0 * y1) - (x1 * y0)
result_x += (x0 + x1) * cross
result_y += (y0 + y1) * cross
result_x /= (area * 6.0)
result_y /= (area * 6.0)
return (result_x, result_y)
import shapefile
sf = shapefile.Reader("cb_2013_us_zcta510_500k.shp")
shapes = sf.iterShapes()
records = sf.iterRecords()
for sh, rec in zip(shapes, records):
print("zipcode: %s" % rec[0])
print(centroid_of_polygon(sh.points))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment