Census maps

using library(tidycensus), made a couple census maps for variables such as: Median household income in the past 12 months (B19013_001), Estimate!!Total:!!Black or African American alone (B02001_003), Estimate!!Total Hispanic or Latino Origin by Race (B03002_001), Estimate!!Total:!!Male: (B01001_002), Estimate!!Total:!!female: (B01001_026).

library(tidycensus)
library(sf)
library(tigris)
library(showtext)
library(sysfonts)
options(tigris_use_cache = TRUE)
# set up with API
census_api_key("XXXXXXXXXXXXXXXXXXXXXXXXX", install = TRUE)
readRenviron("~/.Renviron")


# Based off Sciencbase Map (https://www.sciencebase.gov/catalog/item/imap/63f79b64d34e4f7eda456552), states included in study are: 
# Washington, Oregon, California, Idaho, Nevada, Utah, Arizona, Montana, Wyoming, Colorado, New Mexico, North Dakota, South Dakota, Nebraska, Kansas, Oklahoma, Texas, Minnesota, Iowa, Missouri, Arkansas,  Louisiana

# Parts of Western Wisconsin, Illinois, Tennesse, Mississippi 
# Excluding these states for now until I get more info 

# Take a peak at the available census variables
view_vars <- tidycensus::load_variables(year = 2022, dataset = "acs1") # American Community Survey
view_vars
# nrow(view_vars) #36607

states <- c('Washington', 'Oregon', 'California', 'Idaho', 'Nevada', 'Utah', 'Arizona', 'Montana', 'Wyoming', 'Colorado', 'New Mexico', 'North Dakota', 'South Dakota', 'Nebraska', 'Kansas', 'Oklahoma', 'Texas', 'Minnesota', 'Iowa', 'Missouri', 'Arkansas', 'Louisiana')

proj <- '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs'

counties_outline_col = "grey80"
conus_outline_col = 'grey50'
dpi = 300
bg_col = "white"

font_legend <- 'Source Sans Pro'
font_add_google(font_legend)
showtext_opts(dpi = 300, regular.wt = 200, bold.wt = 700)
showtext_auto(enable = TRUE)

conus_sf <- tigris::states(cb = TRUE) %>%
  st_transform(proj) %>%
  mutate(group = case_when(
    STUSPS %in% c(state.abb[!state.abb %in% c('AK', 'HI')], 'DC') ~ 'CONUS',
    STUSPS %in% c('GU', 'MP') ~ 'GU_MP',
    STUSPS %in% c('PR', 'VI') ~ 'PR_VI',
    TRUE ~ STUSPS
  )) %>%
  filter(group %in% c('CONUS')) |> 
  rmapshaper::ms_simplify(keep = 0.2) |> 
  filter(NAME %in% states)

# Get census data of interest
get_census <- function(geography, variable, year) {
  # Your code for fetching ACS data and processing it
  df <- get_acs(
    geography = geography,
    variable = variable,
    state = states,
    year = year,
    geometry = TRUE
  ) |>
    janitor::clean_names() |>
    st_transform(proj)
  return(df)
}

med_income <- get_census(geography = "county", variable = "B19013_001", year = 2022) 
tot_black <- get_census(geography = "county", variable = "B02001_003", year = 2022) 
tot_latino<- get_census(geography = "county", variable = "B03002_001", year = 2022) 
tot_male <- get_census(geography = "county", variable = "B01001_002", year = 2022)  # Estimate!!Total:!!Male:
tot_female <- get_census(geography = "county", variable = "B01001_026", year = 2022)  # Estimate!!Total:!!Female:

#  Median household income in the past 12 months (B19013)
med_income_map <- med_income |> 
  ggplot(aes(fill = estimate)) + 
  geom_sf(color = counties_outline_col,
          linewidth = 0.05) +
  geom_sf(data = conus_sf,
          fill = NA,
          color = conus_outline_col,
          linewidth = 0.2,
          linetype = "solid") + 
  theme_void() +
  scale_fill_distiller(
    palette = "OrRd", 
    direction = 1,
    name = "Median household income\n2022",
    limits = c(0, max(med_income$estimate)),
    labels = scales::comma) +
  theme(text = element_text(family = font_legend, size = 10),
        legend.margin = margin(r = 10))

ggsave("med_income_census_2022.png", med_income_map, width = 5, height = 4, dpi = dpi, bg =  bg_col)