US Census incorporated places in PostGIS

census_places.py

import urllib2
import sys, os
from StringIO import StringIO
from zipfile import ZipFile
import psycopg2
import subprocess

# Connect to Postgres
pg_conn = psycopg2.connect(dbname="database", user="you", host="localhost", port=5432)
pg_cur = pg_conn.cursor()

pg_cur.execute("""
  DROP TABLE IF EXISTS places;
  CREATE TABLE places (
    statefp character varying(2),
    placefp character varying(5),
    placens character varying(8),
    affgeoid character varying(16),
    geoid character varying(7),
    name character varying(100),
    lsad character varying(2),
    aland double precision,
    awater double precision,
    geom geometry
);
""")
pg_conn.commit()

pg_cur.execute("select substring(fips, char_length(fips) - 1) AS codes FROM us_states")
codes = pg_cur.fetchall()

for code in codes:
  print code[0]
  request = urllib2.urlopen("http://www2.census.gov/geo/tiger/GENZ2013/cb_2013_" + str(code[0]) + "_place_500k.zip")

  with ZipFile(StringIO(request.read())) as zf:
    zf.extractall("shapefiles")

  subprocess.call("ogr2ogr shapefiles/" + str(code[0]) + "_reprojected.shp -t_srs 'EPSG:4326' shapefiles/cb_2013_" + str(code[0]) + "_place_500k.shp", shell=True)

  subprocess.call("shp2pgsql -s 4326 -a shapefiles/" + str(code[0]) + "_reprojected.shp places | psql -U john geomacro", shell=True)

subprocess.call("rm shapefiles/*", shell=True)

SELECT *
FROM places
ORDER BY geom <-> 'SRID=4326;POINT(-89.4 43.07)' limit 10;