Lab8

This map displays the population changes for Rutherford, Williamson,and Davidson between the years 2018-2023 in those TN counties.

Grey and Blue map display the current estimates and the changes.

Code:

# Installing required packages

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)

# Transmitting your Census API key 

census_api_key("5ef79a6f4380ae0dfefc4c256e98d623316f159a")

# Getting ACS variable info
# Change "2023" to your preferred year.

DetailedTables <- load_variables(2023, "acs5", cache = TRUE)
SubjectTables <- load_variables(2023, "acs5/subject", cache = TRUE)
ProfileTables <- load_variables(2023, "acs5/profile", cache = TRUE)

# Using variable DP05_0001 - Total population

# Map data for Period 1

MapData1 <- get_acs(geography = "county subdivision",
                    state = "TN",
                    county = c("Rutherford","Davidson","Williamson"),
                    variables = c(Variable1 = "DP05_0001"),
                    year = 2018,
                    survey = "acs5",
                    output = "wide",
                    geometry = TRUE)

# Map data for Period 2

MapData2 <- get_acs(geography = "county subdivision",
                    state = "TN",
                    county = c("Rutherford","Davidson","Williamson"),
                    variables = c(Variable2 = "DP05_0001"),
                    year = 2023,
                    survey = "acs5",
                    output = "wide",
                    geometry = TRUE)

mapviewOptions(basemaps.color.shuffle = FALSE)

# Calculating change

MapData1_nogeo <- st_drop_geometry(MapData1)
ChangeData <- merge(MapData2, MapData1_nogeo, by = "GEOID")
ChangeData <- ChangeData %>% 
  rename(Area = NAME.x,
         Before = Variable1E,
         Before_M = Variable1M,
         Now = Variable2E,
         Now_M = Variable2M) %>%
  mutate(Change = Now - Before,
         Sig = significance(Before,
                            Now,
                            Before_M,
                            Now_M,
                            clevel = 0.90),
         Sig = case_when(Sig == "TRUE" ~ "Significant",
                         Sig == "FALSE" ~"Nonsignificant",
                         TRUE ~ "Error")) %>% 
  select(Area, Before, Before_M, Now,
         Now_M, Change, Sig, geometry)

# Formatting the blue map

CurrentMap <- mapview(
  MapData2,
  zcol = "Variable2E",
  col.regions = brewer.pal(9, "Blues"),
  layer.name = "Current estimate",
  popup = FALSE)

# Formatting the 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", "Sig")))

# Putting the blue and gray maps on top of each other

CurrentMap | ChangeMap