for_hilary.R

library(grid)
library(maps)
require(ggmap)
library(vegan)
require(mapproj)
library(ggplot2)
library(maptools)
library(pheatmap)
library(gridExtra)

#######################################################################################################################
#######################################################################################################################
#
#   GENERATING THE MAP
#
#######################################################################################################################
#######################################################################################################################

# adapted from http://loloflargenumbers.com/blog/?p=206
us50_shp <- readShapePoly("~/states.shp")
us50_df <- as.data.frame(us50_shp)
us50_points <- sp2tmap(us50_shp)
names(us50_points) <- c("id", "x", "y")
us50 <- merge(x = us50_df, y = us50_points, by.x = "DRAWSEQ", by.y = "id")
cont_us <- us50[us50$STATE_ABBR != "HI" & us50$STATE_ABBR != "AK", ]
ak <- us50[us50$STATE_ABBR == "AK", ]
hi <- us50[us50$STATE_ABBR == "HI", ]

centerState <- function(.df) {
	.df$x <- .df$x - (diff(range(.df$x, na.rm = T))/2 + min(.df$x, na.rm = T))
	.df$y <- .df$y - (diff(range(.df$y, na.rm = T))/2 + min(.df$y, na.rm = T))
	return(.df)
}


scaleState <- function(.df, scale_matrix, scale_factor, x_shift, y_shift) {
	.df <- centerState(.df)
	coords <- t(cbind(.df$x, .df$y))
	scaled_coord <- t(scale_factor*scale_matrix %*% coords)
	.df$x <- scaled_coord[,1] + x_shift
	.df$y <- scaled_coord[,2] + y_shift
	return(.df)
}


update_ak_hi_lat_lon <- function(m){
	# this function takes in a dataframe with lat and lon fields, scales them according to the
	# special placement of AK and HI states on the map.
	for(sample in row.names(m[m$state == "HI", ])){
		hi_point <- data.frame(y=m[row.names(m) == sample, ]$lat,
				x=m[row.names(m) == sample, ]$lon,
				DRAWSEQ = NA, STATE_NAME = "Hawaii",
				STATE_FIPS = 15, SUB_REGION = "point", STATE_ABBR = "HI")
		sp <- scaleState(rbind(hi, hi_point), scale_mat, 1.5, x_shift = -107, y_shift = 25)
		m[row.names(m) == sample, ]$lat <- sp[sp$SUB_REGION == 'point', ]$y
		m[row.names(m) == sample, ]$lon <- sp[sp$SUB_REGION == 'point', ]$x
	}
	
	for(sample in row.names(m[m$state == "AK", ])){
		ak_point <- data.frame(y=m[row.names(m) == sample, ]$lat,
				x=m[row.names(m) == sample, ]$lon,
				DRAWSEQ = NA, STATE_NAME = "Alaska",
				STATE_FIPS = 15, SUB_REGION = "point", STATE_ABBR = "AK")
		sp <- scaleState(rbind(ak, ak_point), scale_mat, 0.4, x_shift = -120, y_shift = 25)
		m[row.names(m) == sample, ]$lat <- sp[sp$SUB_REGION == 'point', ]$y
		m[row.names(m) == sample, ]$lon <- sp[sp$SUB_REGION == 'point', ]$x
	}
	
	return(m)
	
}


scale_mat <- matrix(c(1,0,0,1.25), ncol = 2, byrow = T)
ak_scale <- scaleState(ak, scale_mat, 0.4, x_shift = -120, y_shift = 25)
hi_scale <- scaleState(hi, scale_mat, 1.5, x_shift = -107, y_shift = 25)
all_us <- rbind(cont_us, ak_scale, hi_scale)
all_us$Year.Reported <- NA

WNSData <- read.table(file = '~/WNSData.txt', header = TRUE , sep = "\t" , quote = "")

for(ST in levels(all_us$STATE_ABBR)){
	all_us[all_us$STATE_ABBR == ST, ]$Year.Reported <- WNSData[WNSData$State == ST, ]$Year.Reported
}

pp <- function(d){
	p <- ggplot() + geom_polygon(data=d, aes(x=x, y=y, group = DRAWSEQ, fill=Year.Reported), alpha=0.8, colour='#000000', size=0.2)
	p <- p + theme_bw()
	p <- p + theme(line = element_blank(), line = element_blank(), title = element_blank(), legend.position = "right")

	print(p) 
}

pdf("~/mapd.pdf", width=14, height=8)
pp(all_us)
dev.off()