using tiger census data at the block level to find the population of a circular area

plot_shp.py

"""plot_shp.py

Plots a shapefile with fiona, decartes and matplotlib.
"""

import descartes
import fiona
from matplotlib import pyplot


infile = 'nearby_blocks/nearby_blocks.shp'
savepath = 'out.png'
figure = pyplot.figure(dpi=90)
figure.set_size_inches(9, 9)
axes = figure.add_subplot(1, 1, 1)
blue = '#6699cc'

# Add a patch for each feature in the shapefile
with fiona.collection(infile) as inshp:
  for feature in inshp:
    axes.add_patch(descartes.PolygonPatch(feature['geometry'], fc=blue,
                                          ec=blue, alpha=0.5))
  bounds = inshp.bounds
  axes.set_xlim(bounds[0], bounds[2])
  axes.set_ylim(bounds[1], bounds[3])
  figure.savefig(savepath)

population script output

census blocks in the 2800m search radius: 1061
total number of people in the search radius: 70912
all this took 43.8 seconds

population_within_radius.py

"""population_within_radius.py

 * counts the number of census blocks within a search radius
 * and counts the total number of people living in these blocks
 * saves a new shapefile with just those blocks inside the radius
"""

import math
import time

import fiona

# Assumes you've downloaded the appropriate state's shapefile via
# http://www2.census.gov/geo/tiger/TIGER2010BLKPOPHU
infile = 'tabblock2010_06_pophu/tabblock2010_06_pophu.shp'
outfile = 'nearby_blocks/nearby_blocks.shp'

# We'll use the odd (lon, lat) convention set by the tiger file.
search_center = (-122.081945, 37.397034)  # shoreline and central
search_radius = 2800  # meters


def haversine(point1, point2):
  """Computes distance between two (lon, lat) points in meters."""
  lon1, lat1 = point1
  lon2, lat2 = point2
  lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
  delta_lon = lon2 - lon1
  delta_lat = lat2 - lat1
  a = (math.sin(delta_lat / 2)**2 + math.cos(lat1) * math.cos(lat2) *
       math.sin(delta_lon / 2)**2)
  c = 2 * math.asin(math.sqrt(a))
  earth_radius = 6367000  # meters
  return earth_radius * c


blocks_in_radius = 0
population_in_radius = 0
start = time.time()

with fiona.open(infile) as source, fiona.open(outfile, 'w', crs=source.crs,
                                              driver=source.driver,
                                              schema=source.schema) as sink:
  # Each record in the source is a census block.
  # We'll iterate through all of them..
  for block in source:
    # Test to see if a random point from the block's perimeter is within the
    # search radius.
    #
    # MultiPolygon geometries have coordinates that are lists of lists of
    # tuples while Polygon geometries have coordinates that are simply lists of
    # tuples.
    if block['geometry']['type'] == 'Polygon':
      sample_block_coordinate = block['geometry']['coordinates'][0][0]
    if block['geometry']['type'] == 'MultiPolygon':
      sample_block_coordinate = block['geometry']['coordinates'][0][0][0]

    distance = haversine(sample_block_coordinate, search_center)
    if distance <= search_radius:
      blocks_in_radius += 1
      population_in_radius += block['properties']['POP10']
      sink.write(block)

print 'census blocks in the %dm search radius: %d' % (search_radius,
                                                      blocks_in_radius)
print 'total number of people in the search radius: %d' % population_in_radius
print 'all this took %0.1f seconds' % (time.time() - start)