library(ggplot2)
library(maps)
library(mapproj)

# Get counties and matching (relevant) series IDs
# (area type == "F" selects county-level data,
#  measure code == "03" selects unemployment rate data,
#  seasonal == "U" selects non-seasonally-adjusted data)
areas <- read.delim("areas", header = F, stringsAsFactors = F, row.names = NULL, strip.white = T)
names(areas) <- c("area_type_code", "area_code", "area_text", "display_level", "selectable", "sort_sequence", NULL)
series_ids <- read.delim("series_ids", stringsAsFactors = F, row.names = NULL, strip.white = T)
aux_info <- merge(areas, series_ids, by = c("area_code","area_type_code"))
aux_info <- subset(aux_info,
                   area_type_code == "F" & measure_code == 3 & seasonal == "U",
                   select = c(series_id, area_code, area_text, area_type_code, measure_code, seasonal))

# Read in BLS data and reduce to county-level,
# non-seasonally-adjusted, annual average
# unemployment rate data.
unemp_05_09 <- read.delim("unemployment_05-09", stringsAsFactors = F, row.names = NULL, strip.white = T)
unemp_10_11 <- read.delim("unemployment_10-14", stringsAsFactors = F, row.names = NULL, strip.white = T)
unemp_full <- rbind(unemp_05_09, unemp_10_11)
unemp_full <- subset(unemp_full, period == "M13", select = c(series_id, year, value, period))
unemp_full <- merge(unemp_full, aux_info, by = "series_id")
unemp_full <- transform(unemp_full, value = as.numeric(value))

# Parse out county names and state abbreviations.
unemp_full$county <- tolower(gsub(" County, [A-Z]{2}", "", unemp_full$area_text))
unemp_full$county <- gsub("^(.*) (city|parish), ..$","\\1", unemp_full$county)
unemp_full$state <- gsub("^.*([A-Z]{2}).*$", "\\1", unemp_full$area_text)

# Bucket rate data.
rate_breaks = c(seq(0, 10, by = 2), 15, 20, 40)
unemp_full$value_discrete <- cut(unemp_full$value, breaks = rate_breaks)

# Get state border and county region data.
states <- map_data("state")
counties <- map_data("county")
names(counties) <- c("long", "lat", "group", "order", "state_name", "county")
counties$state <- state.abb[match(counties$state_name, tolower(state.name))]
counties <- subset(counties, select = -c(state_name))

# Break data out into years.
# For each, merge in county region data, then plot:
# fill regions, then overlay county and state borders.
for (yr in c("2008", "2009", "2010", "2011")) {
  unemp <- subset(unemp_full, year == yr)
  unemp <- merge(counties, unemp, by = c("state", "county"))
  unemp <- unemp[order(unemp$order), ]
  choropleth <- ggplot(unemp, aes(long, lat, group = group)) +
    coord_map(project="globular") +
    geom_polygon(aes(fill = value_discrete)) +
    scale_fill_brewer(palette = "OrRd", breaks = rate_breaks) +
    geom_path(data = counties, color = "#222222", size = 0.2, alpha = 0.5) +
    geom_path(data = states, color = "#555555", size = 0.3) +
    labs(title = yr) +
    theme(panel.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          plot.background = element_blank(),
          axis.ticks = element_blank(),
          axis.text.x = element_blank(), axis.text.y = element_blank(),
          axis.title.x = element_blank(), axis.title.y = element_blank())
  ggsave(filename = paste(c(yr, ".png"), collapse = ""), plot = choropleth)
}