bay_creek_data/bay_creek_data.R

# ---- load libraries
library(tidyverse)
library(sf)
library(ggmap)
library(scales)
library(ggrepel)

setwd("~/Documents/Bay_Creek/bay_creek_data")

# ---- load data
block_data_2022 <- sf::read_sf("data/nip_bg_22/nip_bg_22.shp")
metadata_2022 <- read_csv("data/nip_bg_22/nip_metadata_22.csv")
extent <- st_bbox(block_data_2022)

# ---- define areas of interest
block_groups <- data.frame(name = c("Bay Creek 1", "Bay Creek 2"), geo_id = c("550250013001", "550250013002"))

centroids <- block_data_2022 %>% 
  left_join(block_groups, by = "geo_id") %>% 
  filter(geo_id %in% block_groups$geo_id) %>% 
  st_centroid() %>% 
  pull(geometry) %>%
  transpose()

block_groups["lon"] <- unlist(centroids[[1]])
block_groups["lat"] <- unlist(centroids[[2]])

# ---- download basemap
zoom_level <- 12
if(file.exists(paste0("data/basemap_cache/basemap_", zoom_level, ".RData"))){
  load(file = paste0("data/basemap_cache/basemap_", zoom_level, ".RData"))
} else {
  register_stadiamaps(substr(read_file("data/stadia_api_key.txt"), 1, 36), 
                      write = FALSE)
  basemap <- ggmap::get_stadiamap(bbox = c(left = as.double(extent[1]), 
                                           bottom = as.double(extent[2]), 
                                           right = as.double(extent[3]), 
                                           top = as.double(extent[4])),
                                  zoom = zoom_level,
                                  maptype = "alidade_smooth",
                                  color = "bw",
                                  force = TRUE)
  
  save(basemap_raster, file = paste0("data/basemap_cache/basemap_", zoom_level, ".RData"))
}

# ---- plot figures
# --- plot maps
ggmap(basemap) +
  geom_sf(data = block_data_2022,
               aes(fill = medhhinc),
          inherit.aes = FALSE,
          alpha = 0.6) +
  geom_label_repel(data = block_groups,
    aes(label = name, 
        y = lat,
        x = lon),
    min.segment.length = 0,
    nudge_y = -0.03) +
  # geom_sf_label(data = block_data_2022 %>% left_join(block_groups, by = "geo_id"),
  #               aes(label = name),
  #               inherit.aes = FALSE,
  #               nudge_x = 1,
  #               size = 2) +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank()) +
  scale_fill_continuous(label = scales::label_dollar(),
                        type = "viridis") +
  labs(title = "Median Income",
       fill = NULL)
ggsave(file = "figures/median_income_map.png",
       device = "png",
       width = 11,
       height = 8.5,
       units = "in")

# ---- plot graphs
ggplot(data = block_data_2022 %>% mutate(baycreek = geo_id %in% block_groups$geo_id) %>% left_join(block_groups, by = "geo_id") %>% arrange(medhhinc)) +
  geom_hline(data = block_data_2022 %>% filter(geo_id == "Madison"),
             aes(yintercept = medhhinc),
             linetype = "dashed") +
  geom_label(data = block_data_2022 %>% filter(geo_id == "Madison"),
             aes(y = medhhinc,
                 x = 20),
             label = "Madison median") +
  geom_col(aes(x = reorder(geo_id, medhhinc, sum),
               y = medhhinc,
               fill = baycreek)) +
  geom_label(aes(x = reorder(geo_id, medhhinc, sum),
                 y = medhhinc + 10000,
                 label = name)) +
  scale_x_discrete(labels = NULL, breaks = NULL) +
  scale_y_continuous(label = scales::label_dollar()) +
  scale_fill_discrete(guide="none") +
  theme(axis.text.x=element_blank(),
        axis.title.x=element_blank()) +
  labs(title = "Median Income by Block Group",
       x = NULL,
       y = NULL)
ggsave(file = "figures/median_income.png",
       device = "png",
       width = 11,
       height = 8.5,
       units = "in")