# Packages need to be installed (you only need to do this once)

install.packages("maps")
install.packages("mapproj")

# Function for mapping the 48 states

statemaps <- function (a, grayscale=FALSE, ...){
  if (length(a)==51){
    no.dc <- c(1:8,10:51)
    a <- a[no.dc]
  }
  if (length(a)==50){
    lower48 <- state.abb!="AK" & state.abb!="HI"
    a <- a[lower48]
  }
  else if (length(a)!=48) stop ("wrong number of states")
  mapping <- list (1,2,3,4,5,6,7,9,10,11,12,13,14,15,16,17,18,19,20:22,23:24,25,26,27,28,29,30,31,32,33,34:37,38:40,41,42,43,44,45,46,47,48,49,50,51,52,53:55,56:60,61,62,63)
  # for (i in 1:length(mapping)) print(regions[mapping[[i]]])
  a.long <- rep (NA, 63)
  projection <- "bonne"
  for (i in 1:48){
    a.long[mapping[[i]]] <- a[i]
  }
  if (grayscale){
    a.long.scaled <- .95*(a.long-min(a,na.rm=TRUE))/(max(a,na.rm=TRUE)-min(a,na.rm=TRUE))
    shades <- a.long.scaled
    not.dc <- !is.na(a.long.scaled)
    shades[not.dc] <- gray (shades[not.dc])
    map('state', proj=projection, param=25, lty=0, ...)
    m <- map('state', proj=projection, param=25, fill=TRUE, plot=FALSE)
    polygon(m$x,m$y, col=shades, lwd=0.5, border="gray30")
  }
  else {
    map('state', proj=projection, param=25, lty=0, ...)
    m <- map('state', proj=projection, param=25, fill=TRUE, plot=FALSE)
    polygon(m$x,m$y, col=a.long, lwd=0.5, border="gray30")
  }
}

# Functions that make a pretty color scheme

blrd <- function(n) {
  if (n %% 2 == 1){
    n2 <- (n + 1)/2
    return(c(hsv(4/6,1,0.6,alpha=seq(1,0,length.out=n2))[1:(n2-1)], "white", hsv(1,1,0.6,alpha=seq(0,1,length.out=n2))[2:n2]))
  }
  else if (n %% 2 == 0) {
    n2 <- n/2
    return(c(hsv(4/6,1,0.6,alpha=seq(1,0,length.out=(n2+1)))[1:n2], hsv(1,1,0.6,alpha=seq(0,1,length.out=(n2+1)))[2:(n2+1)]))
  }
}

redblue <- function (a, lo=min(a), hi=max(a), n.shades=1001){
  a <- pmin (hi, pmax (lo, a)) 
  return (blrd(n.shades)[floor (1 + (n.shades-1)*(a-lo)/(hi-lo))])
}

# Data from 2004 and 2008 elections

elections <- read.table("http://www.stat.columbia.edu/~gelman/surveys.course/pres20042008.txt")
no_dc <- rownames(elections) != "DC"

# Now let's make some maps!

png ("2008_2004_residuals.png", height=300, width=400)
par (mar=c(2,0,0,0))
M1 <- lm (obama08 ~ kerry04, data=elections, subset=no_dc)
statemaps (resid(M1), grayscale=TRUE)
mtext ("Residuals from predicting Obama08 from Kerry04", line=-1)
mtext ("Light colors:  Obama did relatively well in 2008,\nDark colors:  McCain did relatively well in 2008", side=1, line=0, cex=.8)
dev.off()

png ("2008_2004_difference.png", height=300, width=400)
par (mar=c(0,0,0,0))
x <- elections[no_dc,"kerry04"]
y <- elections[no_dc,"obama08"]
statemaps (redblue(x-y, -.1, .1))
mtext ("Obama did better than Kerry in nearly every state", line=-2)
mtext ("Blue:  Obama did better in 2008 than Kerry did in 2004,\nRed:  Obama did worse than Kerry,\nLight colors:  Little or no change", side=1, line=-1, cex=.8)
dev.off()