generate_us_map_for_sewage_study.R

library(ggplot2)
library(maps)
require(ggmap)
library(vegan)
require(mapproj)
library(maptools)

#######################################################################################################################
#######################################################################################################################
#
#   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)
}

centerState(ak)$y
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)
}

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)

# file to map states into climate regions
# http://www.ncdc.noaa.gov/monitoring-references/maps/us-climate-regions.php#references
# Thomas R. Karl and Walter James Koss, 1984: "Regional and National Monthly, Seasonal, and Annual Temperature Weighted by Area, 1895-1983." Historical Climatology Series 4-3, National Climatic Data Center, Asheville, NC, 38 pp.
state_region <- read.table(file = '~/00_state_region_data.txt', header = TRUE , sep = "\t" , quote = "")

# replace SUB_REGION with actual climate regions in the US
levels(all_us$SUB_REGION) <- c(levels(all_us$SUB_REGION), levels(state_region$region))
for(ST in levels(all_us$STATE_ABBR)){
	all_us[all_us$STATE_ABBR == ST, ]$SUB_REGION <- state_region[state_region$state == ST, ]$region
}

# metadata file for sewage dataset
metadata <- read.table(file = '~/00_metadata.txt', header = TRUE , sep = "\t" , quote = "")

for(sample in metadata[metadata$state == "HI", ]$samples){
	hi_point <- data.frame(y=metadata[metadata$samples == sample, ]$lat, x=metadata[metadata$samples == 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)
	metadata[metadata$samples == sample, ]$lat <- sp[sp$SUB_REGION == 'point', ]$y
	metadata[metadata$samples == sample, ]$lon <- sp[sp$SUB_REGION == 'point', ]$x
}

for(sample in metadata[metadata$state == "AK", ]$samples){
	ak_point <- data.frame(y=metadata[metadata$samples == sample, ]$lat,
							x=metadata[metadata$samples == 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)
	metadata[metadata$samples == sample, ]$lat <- sp[sp$SUB_REGION == 'point', ]$y
	metadata[metadata$samples == sample, ]$lon <- sp[sp$SUB_REGION == 'point', ]$x
}


# unique city names to avoid overprinting:
cdf <- metadata[!duplicated(metadata[,c('city')]),]
# unique latitudes to avoid overprinting:
pdp <- metadata[!duplicated(metadata[,c('lat')]),]

pp <- function(d){
	p <- ggplot() + geom_polygon(data=d, aes(x=x, y=y, group = DRAWSEQ, fill=SUB_REGION), alpha=0.2, colour='#000000', size=0.2)
	p <- p + geom_text(data=cdf, aes(x=lon, y=lat + 0.8, label=city), colour="#000000", alpha=0.4)
	p <- p + geom_point(data=pdp, aes(x=lon, y=lat), size = 5, colour="#FF0000", alpha=0.15)
	p <- p + geom_point(data=pdp, aes(x=lon, y=lat), size = 3, colour="#FF0000", alpha=0.75)
	p <- p + geom_point(data=pdp, aes(x=lon, y=lat), colour="#000000", size = 1, alpha=0.75)
	p <- p + theme_bw()
	p <- p + theme(line = element_blank(), text = element_blank(), line = element_blank(), title = element_blank(), legend.position = "none")
	print(p) 
}

# yay for the map!
pdf("~/map.pdf", width=12, height=8)
pp(all_us)
dev.off()