Amanda Cox's R sketches for NYT drought map

drought.R

# code by Amanda Cox, shared at NICAR13 (http://d.pr/f/D2yw)
# http://www.nytimes.com/interactive/2012/07/20/us/drought-footprint.html


#install.packages("gdata")
library(gdata)
library(maptools)
library(RColorBrewer)

### set our working directory
#setwd("/Users/192363/Desktop/droughtAC")
setwd("C://training//NYT-R//droughtAC")


### read some data
data <- read.xls("drought.xls")
map <- readShapePoly("drought.shp")

### look at what we got.
plot(map)
plot(map, col = "blue", border = "white")

names(data)
nrow(data)
head(data)
data$Y2012
hist(data$Y2012)
?hist
hist(data$Y2012, breaks = 50, col = "lightyellow")
abline(v = -1.25, col = "red")

### super important. make our data in the same order as the shapes in the map.
data <- data[match(map$AREA, data$AREA), ]

### add a color field to our data
data$color <- ifelse(data$Y2012 < -1.25, "red", "lightgrey")
data$color
table(data$color)

## make a simple map
plot(map, col = data$color)

#### learn about for loops
for(i in 1:10){
  print(i)
}

### make a bunch of maps in a for loop
par(mfrow = c(4,4))
for(i in 2000:2012){
	var <- paste0("Y", i)
	color <- ifelse(data[ ,var] <= -1.25, "red", "lightgrey")
	plot(map, col = color, border = F)
	title(var)
}


### let's make our plotting a function
drawMaps <- function(years, threshold = -1.25){
	for(i in years){
		var <- paste0("Y", i)
		color <- ifelse(data[ ,var] <= threshold, "red", "lightgrey")
		plot(map, col = color, border = F)
		title(var)
	}	
}

### make a bunch of maps
pdf(file = "droughtmaps.pdf", ,width=8,height=11)
	par(mfrow=c(12,10), mar=c(0,0,0,0))
	drawMaps(1900:2012)
dev.off()


### let's get fancier with our colors
install.packages("RColorBrewer")
library(RColorBrewer)
help.search("RColorBrewer")
display.brewer.all(n=10, exact.n=FALSE)

mycolors <- brewer.pal(5, "YlOrRd")
mycolors <- rev(mycolors)
plot(map, col = mycolors[cut(data$Y2012, breaks = 5)])


### let's make line charts instead
plot(data$Y2012, type = "l")

years <- 1895:2012

plot(t(data[1, paste0("Y", years)]), type = "l")
lines(t(data[2, paste0("Y", years)]), col = "red")

#### lets make a fancier line chart
toPlot <- data[, paste0("Y",years)]
plot(0, type = "n", ylim = range(toPlot), xlim=c(1895,2012), xlab = "", ylab = "")
apply(toPlot, 1, function(x) lines(years, x, col = adjustcolor("lightgrey", alpha = .05)))
worst <- which.min(data$Y2012)
lines(years, t(data[worst, paste0("Y", years)]), type="l", col = "red")
text(2012, data$Y2012[worst], paste(data$NAME, data$ST.C.2)[worst], adj = 1)

### let's compute some percentages by state
data$area <- map$AREA
data$drought <- map$AREA * as.numeric(data$Y2012<c(-1.25))
bystate <- aggregate(data[,c("area", "drought")], list(state = data$ST.C.2), sum)
bystate$pctdrought <- 100*bystate$drought/bystate$area
bystate[order(bystate$pctdrought),]

### let's do some weird clustering just for fun
z <- cmdscale(dist(toPlot))
plot(z, cex = .25)
text(z[,1], z[,2], paste(casefold(as.character(data$NAME)), data$ST.C.2), cex  = .5)

### hardcore friends only: let's reproject our map
library(rgdal)
aea.proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m"  # comes from the internet (http://lists.maptools.org/pipermail/proj/2003-March/000665.html)
ny.proj <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +units=m"
proj4string(map) = CRS(aea.proj)
map2 <- spTransform(map, CRS(ny.proj))