ashaw · March 3, 2013 16:16
diff --git a/drought.R b/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))
	# 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))