Lab 8 Map

Population from 2018-2023

This interactive map shows the population in Rutherford, Davidson, and Williamson County. Simply use the slider to see the change from 2018-2023:

Code:

# Load required libraries
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("tidycensus")) install.packages("tidycensus")
if (!require("sf")) install.packages("sf")
if (!require("mapview")) install.packages("mapview")
if (!require("leaflet.extras2")) install.packages("leaflet.extras2")
if (!require("scales")) install.packages("scales")
if (!require("RColorBrewer")) install.packages("RColorBrewer")

library(tidycensus)
library(tidyverse)
library(sf)
library(mapview)
library(leaflet.extras2)
library(scales)
library(leaflet)
library(leafpop)
library(RColorBrewer)

options(tigris_use_cache = TRUE)

# Census API key (consider storing this securely instead)
census_api_key("494018a70f114d4f76b10537730ccc9c7dbfe36b")

# Define the ACS variable for total population
VariableList <- c("DP05_0001E") 

# Map data for Period 1 (2018)
MapData1 <- get_acs(
  geography = "county subdivision",
  state = "TN",
  county = c("Rutherford", "Davidson", "Williamson"),
  variables = VariableList,
  year = 2018,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
)

# Map data for Period 2 (2023)
MapData2 <- get_acs(
  geography = "county subdivision",
  state = "TN",
  county = c("Rutherford", "Davidson", "Williamson"),
  variables = VariableList,
  year = 2023,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
)

# Check if all counties are included
unique(MapData2$NAME)  # Should list all 3 counties

# Remove geometry for merging
MapData1_nogeo <- st_drop_geometry(MapData1)  

# Merge data and ensure all areas are included
ChangeData <- merge(MapData2, MapData1_nogeo, by = "GEOID", all.x = TRUE, all.y = TRUE)

# Rename columns for clarity
ChangeData <- ChangeData %>% 
  rename(Area = NAME.x,
         Before = DP05_0001E.y,  # 2018 estimate
         Now = DP05_0001E.x) %>%  # 2023 estimate
  mutate(Change = Now - Before) %>%
  select(Area, Before, Now, Change, geometry)

# Mapping: Current estimate (Blue map)
CurrentMap <- mapview(
  MapData2,
  zcol = "DP05_0001E",
  col.regions = brewer.pal(9, "Blues"),
  layer.name = "Current estimate",
  popup = FALSE
)

# Mapping: Change over time (Gray map)
ChangeMap <- mapview(
  ChangeData,
  zcol = "Change",
  col.regions = brewer.pal(9, "Greys"),
  layer.name = "Change",
  popup = popupTable(
    ChangeData,
    feature.id = FALSE,
    row.numbers = FALSE,
    zcol = c("Area", "Before", "Now", "Change"))
)

# Overlay maps
CurrentMap | ChangeMap