map_response_rate_tract.sh

#!/bin/bash
# Requires, gdal, wget, spatialite, R with (maptools, RColorBrewer, classInt)
# You should be able to put any state/county combination in here and get a map made 
# Replace spaces with an underscore ie:
#state=New_York
#county=New_York
#
#state=ILLINOIS
#county=Cook
#
#state=Indiana
#county=Tippecanoe

state=$1
county=$2

wget -N -nv "http://www.jduck.net.s3.amazonaws.com/files/fips.txt" # Need some data to build ghastly census.gov URLs.
lookup=$(grep -i "${state}|${county}" fips.txt)
if [ -n "${lookup}" ]; then
    stcofips=$(echo $lookup | cut -d"|" -f1) # grab stcofips
    stfips=$(echo ${lookup} | cut -d"|" -f3) # chop out stfips
    cofips=$(echo ${lookup} | cut -d"|" -f2) # chop out cofips 
    STATE=$(echo ${lookup} | cut -d"|" -f4) # state name
    COUNTY=$(echo ${lookup} | cut -d"|" -f5) # county name
else
    echo "Couldn't find that county, did you spell it right?"
    exit
fi

wget -nc "http://www2.census.gov/geo/tiger/TIGER2009/${stfips}_${STATE}/${stcofips}_${COUNTY}/tl_2009_${stcofips}_tract00.zip"

unzip -o tl_2009_${stcofips}_tract00.zip
ogr2ogr -t_srs EPSG:2163 tl_2009_${stcofips}_tract00_aea.shp tl_2009_${stcofips}_tract00.shp 
grep Tract participationrates2010.txt | grep "^${stcofips}" | sed 's/\|\| /\|/g;s/ \|\|/\|/g' > tract_${stcofips}.csv

spatialite << SQL
.mode csv
.separator |
select InitSpatialMetaData();
insert into spatial_ref_sys values (2163,'epsg',2163,'US National Atlas Equal Area','+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs');
create table census_status (fips character(5), tract varchar(50), desc varchar(20), y2000 numeric(5,2), y2010 numeric(5,2));
.import tract_${stcofips}.csv census_status
CREATE virtual table virtual_tract using VirtualShape(tl_2009_${stcofips}_tract00_aea,UTF-8,2163);
create table export as select a.Geometry, a.*, b.* from virtual_tract as a left join (select fips, tract, desc, y2000, y2010 from census_status) as b on b.fips = a.CTIDFP00 where a.TRACTCE00 != "0" and b.y2010 > 0 ;
.dumpshp export Geometry export UTF-8 POLYGON
SQL

export outfile=part_${STATE}_${COUNTY}.pdf
export state=${STATE}
export county=${COUNTY}

R --slave << "EOF"
library(maptools)
library(RColorBrewer)
library(classInt)

outfile <- Sys.getenv(c("outfile"))
state <- Sys.getenv("state")
county <- Sys.getenv("county")

c <- readShapePoly('export.shp')

nclr <- 9
colpal <- brewer.pal(nclr,"Blues")

class <- classIntervals(c$y2010, nclr, style="quantile")
pcol <- findColours(class,colpal)
pdf(outfile, width=17, height=11)
plot(c, col=pcol, border="gray")
title(paste(state, county, sep=" "))
legend("topleft", legend=names(attr(pcol, "table")), fill=attr(pcol, "palette"), cex = 0.75, bty = "n")

dev.off()

EOF