# Use ACS Subject Tables to loand data into  R
# Subject Tables can be found in 
# https://www.census.gov/acs/www/data/data-tables-and-tools/subject-tables/
#Median Income in the Past 12 Months	median income, income	S1903
#  Geography:
# https://walker-data.com/tidycensus/articles/basic-usage.html#geography-in-tidycensus
#  Variables: v17 <- load_variables(2017, "acs5", cache = TRUE)
#  6/05/2020
# Jin Man Lee

# # packages
# install.packages("tidyverse")
# install.packages("tidycensus")
# install.packages("tmap")
# install.packages("tmaptools")

library(tidyverse)
library(tidycensus)
library(ggplot2)
library(tmap)
library(tmaptools)
library(sf)
library(mapview)

#1. Map by census tract 
census_api_key("YOUR KEY GOES HERE")  # Census API key. Obtain one at http://api.census.gov/data/key_signup.html

cook <- get_acs(state = "IL", county = "Cook", geography = "tract", 
                  variables = "B19013_001", geometry = TRUE, year=2018)
head(cook)
mapview(cook, zcol = "estimate", legend = TRUE)

cook %>%
  ggplot(aes(fill = estimate)) + 
  geom_sf(color = "lightblue") + 
  coord_sf() + 
  scale_fill_viridis_c(option = "plasma") 


# Map by census state
state <- get_acs(geography = "state", 
                variables = "B19013_001", geometry = TRUE, year=2018)
head(state)
mapview(state, zcol = "estimate", legend = TRUE)

state %>%
  ggplot(aes(fill = estimate)) + 
  geom_sf() + 
  coord_sf() + 
  scale_fill_viridis_c(option = "magma") 

# Map by County
IL <- get_acs(geography = "county", variables = "B19013_001", state = "IL",year=2018,geometry = TRUE)
mapview(IL, zcol = "estimate", legend = TRUE)

IL %>%
  ggplot(aes(fill = estimate)) + 
  geom_sf() + 
  coord_sf() + 
  scale_fill_viridis_c(option = "plasma") 

IL %>%
  mutate(NAME = gsub(" County, Illinois", "", NAME)) %>%
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) +
  geom_errorbarh(aes(xmin = estimate - moe, xmax = estimate + moe)) +
  geom_point(color = "red", size = 3) +
  labs(title = "Household income by county in Illinois",
       subtitle = "2014-2018 American Community Survey",
       y = "",
       x = "ACS estimate (bars represent margin of error)")


# Map Household by PUMA
cook_puma <- get_acs(geography = "public use microdata area", 
              variables = "B19013_001", state = "IL",year=2018,geometry = TRUE)
cook_puma <- filter(cook_puma,substr(GEOID,1,5)=="17034" |substr(GEOID,1,5)=="17035"  )
mapview(cook_puma, zcol = "estimate", legend = TRUE)

cook_puma %>%
  ggplot(aes(fill = estimate)) + 
  geom_sf() + 
  coord_sf() + 
  scale_fill_viridis_c(option = "magma") 

cook_puma  %>%
  mutate(NAME = gsub(" County, Illinois", "", GEOID)) %>%
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) +
  geom_errorbarh(aes(xmin = estimate - moe, xmax = estimate + moe)) +
  geom_point(color = "red", size = 3) +
  labs(title = "Household income by PUMA in Illinois",
       subtitle = "2014-2018 American Community Survey",
       y = "",
       x = "ACS estimate (bars represent margin of error)")