Map

Population Changes

The interactive map pictured below shows the population changes in Rutherford, Davidson, and Williamson county between 2018 and 2023.

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("2e2a0ce2c78a86d5656ad34469c6cd33e36076b9")

# 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