popBase.R

x = c("ggplot2", "classInt", "dplyr", "rgdal", "maptools", "mapproj", "rgeos", "RCurl")
lapply(x, library, character.only = TRUE)
intervals = function(.df, ...){
	argList = match.call(expand.dots=FALSE)$... 
	for(i in 1:length(argList)){
		colName <- argList[[i]]
		series_colName = eval(substitute(colName), envir=.df, enclos=parent.frame())
		min <- min(series_colName)
		max <- max(series_colName)
		diff <- max - min
		std <- sd(series_colName)
		equal.interval <- seq(min, max, by = diff/6)
		quantile.interval <- quantile(series_colName, probs = seq(0, 1, by = 1/6))
		std.interval <- c(seq(min, max, by = std), max)
		natural.interval <- classIntervals(series_colName, n = 6, style = 'jenks')$brks
		.df$equal <- cut(series_colName, breaks = equal.interval, include.lowest = TRUE)
		names(.df)[names(.df)=="equal"] <- paste(colName,".","equal", sep = '')
		.df$quantile <- cut(series_colName, breaks = quantile.interval, include.lowest = TRUE)
		names(.df)[names(.df)=="quantile"] <- paste(colName,".","quantile", sep = '')
		.df$std <- cut(series_colName, breaks = std.interval, include.lowest =  TRUE)
		names(.df)[names(.df)=="std"] <- paste(colName,".","std", sep = '')
		.df$natural <- cut(series_colName, breaks = natural.interval, include.lowest = TRUE)
		names(.df)[names(.df)=="natural"] <- paste(colName,".","natural", sep = '')
	}
	return(.df)
}
remove.territories = function(.df) {
	subset(.df, 
		.df$id != "AS" &
		.df$id != "MP" &
		.df$id != "GU" & 
		.df$id != "PR" &
		.df$id != "VI" 
	)
}
plain_theme = theme(axis.text=element_blank()) + 
	theme(panel.background = element_blank(), 
		panel.grid = element_blank(), 
		axis.ticks = element_blank())
no_ylab = ylab("") 
no_xlab = xlab("")
pub_date = format(Sys.time(), "%Y-%m-%d")

# See https://rpubs.com/technocrat/thematic-walkthrough for how to create and save cb5
load("data/cb5.Rda")
population<- read.csv(textConnection(getURL('https://www.census.gov/popest/data/national/totals/2014/files/NST_EST2014_ALLDATA.csv')))[-(1:5),]
pop = data.frame(population$STATE,population$POPESTIMATE2014)
names(pop) = c("fips", "pop")
converter = read.csv("https://gist.githubusercontent.com/technocrat/93470bf9abead06ef926/raw/f652f8171374e7808455f42167f5480ea15f7f4e/state_fips_postal.csv", header = FALSE)
converter = rename(converter, state = V1, fips = V2, id = V3)

us = cb5
us_aea = spTransform(us, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
us_aea@data$id = rownames(us_aea@data)
alaska = us_aea[us_aea$STATEFP=="02",]
alaska = elide(alaska, rotate=-50)
alaska = elide(alaska, scale=max(apply(bbox(alaska), 1, diff)) / 2.3)
alaska = elide(alaska, shift=c(-2100000, -2500000))
proj4string(alaska) = proj4string(us_aea)
hawaii = us_aea[us_aea$STATEFP=="15",]
hawaii = elide(hawaii, rotate=-35)
hawaii = elide(hawaii, shift=c(5400000, -1400000))
proj4string(hawaii) = proj4string(us_aea)
us_aea = us_aea[!us_aea$STATEFP %in% c("02", "15"),]
us_aea = rbind(us_aea, alaska, hawaii)
rhode_island = us_aea[us_aea$STATEFP=="44",]
rhode_island = elide(rhode_island, shift=c(0,-200000))
proj4string(rhode_island) = proj4string(us_aea)
delaware = us_aea[us_aea$STATEFP=="10",]
delaware = elide(delaware, shift=c(200000,0))
proj4string(delaware) = proj4string(us_aea)
dc = us_aea[us_aea$STATEFP=="11",]
dc = elide(dc, shift=c(250000,-125000))
proj4string(dc) = proj4string(us_aea)
us_aea = us_aea[!us_aea$STATEFP %in% c("44", "10", "11"),]
us_aea = rbind(us_aea, delaware, rhode_island, dc)
centroids = data.frame(us_aea$STUSPS, coordinates(us_aea))
names(centroids) = c("id", "clong", "clat")
centroids = remove.territories(centroids)
us50 <- fortify(us_aea, region="STUSPS")
us50 = remove.territories(us50)

rownames(centroids) = centroids$id
centroids['LA',]$clong = 729000
centroids['FL',]$clong = 1800000
centroids['DE',]$clong = 2400000
centroids['RI',]$clong = 2400000
centroids['MD',]$clong = 1950000
centroids['MD',]$clat = -360000
centroids['MA',]$clat = 100000
centroids['NJ',]$clong = 2130000
centroids['NJ',]$clat = -255000
  
pop = merge(pop, converter)
pop = intervals(pop, pop)
map = merge(us50,pop) # has to be geo first

b = ggplot(data= map) + 
	geom_map(map=map, aes(x=long, y=lat, map_id=id, group=group), ,fill="white", color="black", size=0.3) + 
	no_ylab + 
	no_xlab + 
	plain_theme

l = geom_text(data=centroids, aes(clong, clat, label = id), color = "black", size=2)


# l is a ggplot ready to add to geom_polygon (must be overplotted)