iPython note:

ZIP-communityWeights

 

This was a favour for the lovely missus – She needed to look at census data , provided by ZIP code, but was working with an organisation that was interested in community areas.

Census data can be obtained in the form of Population/unemployment/Ethnic minority in each zipcode.

To translate this to Community area, we can make an approximation by weighting each zipcode by the portion of it’s Area in each Community area.

This is not a perfect mapping, for example, in zipcode 60101 there may be a part that overlaps Community Area ‘LOOP’.   in this case nobody from 60101 would live in the LOOP, but you would assign some of the unemployed population to that Community area.

Anyway, it’s a first approximation and that’s all we need for now.

Shapefiles for Community Areas are here:

https://data.cityofchicago.org/Facilities-Geographic-Boundaries/Boundaries-Community-Areas-current-/cauq-8yn6

Shapefiles for ZIPs are here:

https://data.cityofchicago.org/Facilities-Geographic-Boundaries/Boundaries-ZIP-Codes/gdcf-axmw

 

Code Highlights osgeo package and Pandas (of course)
from osgeo import ogr
import os
import pandas as pd #My favourite Package!

 

Import a shapefile downloaded from the cityofchicago.org:

driver = ogr.GetDriverByName('ESRI Shapefile')
daShapefile = r"yourCommunitiesFile1.shp"
Communities = driver.Open(daShapefile, 0)
# 0 means read-only. 1 means writeable.

# Check to see if shapefile is found.
if Communities is None:
print 'Could not open %s' % (daShapefile)
else:
print 'Opened %s' % (daShapefile)
com_layer = Communities.GetLayer()
featureCount = com_layer.GetFeatureCount()
print "Number of features in %s: %d" % (os.path.basename(daShapefile),featureCount)

daShapefile2 = r"yourZIPShapeFile1.shp"
ZIP = driver.Open(daShapefile2, 0)
# 0 means read-only. 1 means writeable.
if ZIP is None:
print 'Could not open %s' % (daShapefile2)
else:
print 'Opened %s' % (daShapefile2)
z_layer = ZIP.GetLayer()
featureCount = z_layer.GetFeatureCount()
print "Number of features in %s: %d" % (os.path.basename(daShapefile2),featureCount)

z_layer = ZIP.GetLayer()
c_layer = Communities.GetLayer()

Create a dictionary:

wt = {}
for z_lyr in [z_layer.GetFeature(i) for i in range(z_layer.GetFeatureCount())]:
z = z_lyr.zip
zg = z_lyr.geometry()
wt[z]=[]
    for c_lyr in [c_layer.GetFeature(j) for j in range(c_layer.GetFeatureCount())]:
cg = c_lyr.geometry()
intersection = cg.Intersection(zg)
wt[z].append(intersection.Area()/zg.Area())

Create a list of the communities:

communities = []
for c_lyr in [c_layer.GetFeature(j) for j in range(c_layer.GetFeatureCount())]:
communities.append(c_lyr.community)

 

Create a pandas DataFrame and export to excel:

tbl = pd.DataFrame(wt).T
tbl.columns=communities
tbl.to_excel('zip-community weights.xls')
ZIP-communityWeights