# Toolbox Tuesday 
# This script demonstrates how to extract, process, and map tract-level data for a selected city
# Contact Echo Zheng (zheng@scag.ca.gov) for questions or additional data resources 

#------------------------------------------------------------------------------#
#####                           Analytic Workflow                          #####
#------------------------------------------------------------------------------#
# Analytic goal: Extract, process, and map neighborhood-level (tract) demographic/SES characteristics within a city
# Prerequisites: understand census geography

# (1) Read in a local csv file (a place-to-tract crosswalk)
# (2) Get census data at the tract level (e.g., median household income in 2022 5-year ACS)
# (3) Extract census tracts for a selected city (data wrangling)
# (4) Join city tracts with tract variable(s)
# (5) Map variable across tracts within the city
# (6) Advanced: extract multiple variables 

#-------------------------------------------------------------------------------
#####                         Getting started in R                         #####
#-------------------------------------------------------------------------------

# To get started, install the needed packages - only need to do it once
# Press shift+ctrl+c to comment/un-comment lines
# install.packages("tidycensus")
# install.packages("tidyverse")          # A suite of packages designed for data wrangling
# install.packages("mapview")            # For interactive viewing of spatial data in R
# install.packages("sf")

# Load packages - do it each time you start your analysis
library(tidycensus)
library(tidyverse)
library(mapview)
library(sf)

# Check current working directory
getwd()

# Set your working directory
setwd("C:/Users/zheng/Desktop/ToolboxTuesday2024")
dir.create("output")

# Optional (Recommended): get and activate your Census API key
# visit https://api.census.gov/data/key_signup.html to request a key, then activate the key from the link in your email.
# Set your API key in R
census_api_key("054fae4b3f5a4739d7a9634a2bab59d7b4de32d3", install=TRUE, overwrite = TRUE)

#-------------------------------------------------------------------------------
#####         Read in a local csv file (a place-to-tract crosswalk)        #####
#-------------------------------------------------------------------------------
crosswalk <- read_csv("place-to-tract-crosswalk-SCAG-2020.csv", skip = 1)

#-------------------------------------------------------------------------------
#####                  Get variable(s) from 2022 5-year ACS                #####
#-------------------------------------------------------------------------------
?get_acs()

acsdata <- get_acs(geography = "tract",
                   variables = "B19013_001", 
                   state = "CA",
                   county ="Orange", 
                   year = 2022, 
                   survey = "acs5",
                   geometry = TRUE)  # Set geometry = TRUE to allow for mapping

# Let's map median household income in Orange County!
mapview(acsdata, zcol = "estimate")  

#-------------------------------------------------------------------------------
#####                      Extract tracts within a city                    #####
#-------------------------------------------------------------------------------
# Data wrangling: (1) filter for a subset of data, (2) create new variables, (3) select needed variables
# We use a pipe operator to chain the above three operations 
# Pipe operator: %>% (keyboard shortcuts: Ctrl + Shift + M for Windows or Cmd + Shift + M for Mac)
# Key functions: filter, mutate, rename, select

city_tracts <- crosswalk %>% 
  filter(`Place name` == "Santa Ana city, CA",
         `place-to-tract allocation factor` > 0.01) %>%  # Adds a condition that tract contains at least 1% of the city's population
  mutate(tract1 = substr(`Tract`,1,4),
         tract2 = substr(`Tract`, 6, 7), 
         GEOID = paste0(`County code`, tract1, tract2)) %>% 
  rename(City = `Place name`) %>% 
  select(City, GEOID)

## Notes. Backticks (``) are used to enclose variable names that would otherwise 
##        be invalid due to containing special characters, spaces. 
##        Use "==" to specify values from existing values; Use "=" to create values

#-------------------------------------------------------------------------------
#####                 Join city tracts with tract variables                #####
#-------------------------------------------------------------------------------
# We want to join our ACS data with city tracts 
citydata <- left_join(city_tracts, acsdata, by = "GEOID")

citydata <- citydata %>% 
  rename(`Median HH Income` = estimate) %>% 
  st_as_sf()    # Covert data into an "sf" object (a type of object that contains spatial aspects)


#-------------------------------------------------------------------------------
#####                 Map variable across tracts within city               #####
#-------------------------------------------------------------------------------
# Simple interactive map with mapview()
mapview(citydata, zcol = "Median HH Income")

acsdata <- acsdata %>% rename(`Median HH Income` = estimate)
mapview(acsdata, zcol = "Median HH Income")

#-------------------------------------------------------------------------------
#####                  Export R generated data as shapefile                #####
#-------------------------------------------------------------------------------
st_write(citydata, "output/citydata.shp",  delete_dsn = TRUE)
st_write(acsdata, "output/countydata_OC.shp",  delete_dsn = TRUE)

#-------------------------------------------------------------------------------
#####    ADVANCED 1: extract multiple variables to create local profile    #####
#-------------------------------------------------------------------------------
# Step 1: find the variables in the Census ACS
varlist <- c(population = "B01001_001", 
            med_age = "B01002_001",
            nh_white = "B03002_003",
            foreign = "B05002_013",
            households = "B25003_001",
            renters = "B25003_003",
            workers = "B08301_001",
            wfh = "B08301_021")

acsdata2 <- get_acs(geography = "place",
                    variables = varlist,
                    state = "CA",
                    year = 2022, 
                    survey = "acs5",
                    geometry = FALSE) 


# Step 2: Convert data from "long" form to "wide" form
acsdata2 <- pivot_wider(acsdata2, id_cols = c(GEOID, NAME), names_from=variable, values_from=estimate)

# Step 3: Keep cities of interest and calculate variables needed
citydata2 <- acsdata2 %>% 
  filter(NAME == "Santa Ana city, California" | NAME=="Tustin city, California" | NAME=="Orange city, California") %>% 
  mutate(pct_nhw = nh_white/population * 100,
         pct_fb = foreign/population * 100,
         pct_renters = renters/households * 100,
         pct_wfh = wfh/workers * 100,
         year = "2018-2022 ACS") %>% 
  select(NAME, population, households, workers, starts_with("pct"), year)

# Step 4: export data to excel
write_csv(citydata2, "output/citydata2.csv")

# Optional: extract the same variables from 2017 5-year ACS file
acsdata3 <- get_acs(geography = "place",
                    variables = varlist,
                    state = "CA",
                    year = 2017, 
                    survey = "acs5",
                    geometry = FALSE) 

acsdata3 <- pivot_wider(acsdata3, id_cols = c(GEOID, NAME), names_from=variable, values_from=estimate)

citydata3 <- acsdata3 %>% 
  filter(NAME == "Santa Ana city, California" | NAME=="Tustin city, California" | NAME=="Orange city, California") %>% 
  mutate(pct_nhw = nh_white/population * 100,
         pct_fb = foreign/population * 100,
         pct_renters = renters/households * 100,
         pct_wfh = wfh/workers * 100,
         year = "2013-2017 ACS") %>% 
  select(NAME, population, households, workers, starts_with("pct"), year)

citydata4 <- rbind(citydata2, citydata3)     ## combine "citydata2" and "citydata3" by rows

write_csv(citydata4, "output/citydata4.csv")

#-------------------------------------------------------------------------------
#####                  ADVANCED 2: extract multiple variables              #####
#-------------------------------------------------------------------------------
# Step 1: find the variables in the Census ACS
varlist <- c(workers = "B08301_001", 
             wfh = "B08301_021")

acsdata3 <- get_acs(geography = "tract",  
                   variables = varlist,
                   state = "CA",
                   county =c("Los Angeles", "Orange"),
                   year = 2022, 
                   survey = "acs5",
                   geometry = TRUE)  

# Step 2: Convert data from "long" form to "wide" form
acsdata3 <- acsdata3 %>% as.data.frame()

acsdata3 <- pivot_wider(acsdata3, id_cols = c(GEOID, NAME, geometry), names_from=variable, values_from=estimate)

# Step 3: Calculate percentage of workers working from home
acsdata3 <- acsdata3 %>% 
  mutate(`Percent WFH` = wfh / workers *100) %>% 
  st_as_sf()

#Step 4: Map WFH across Orange County
mapview(acsdata3, zcol = "Percent WFH")


  