Population Statistics Interactive Map: Rutherford, Davidson, and Williamson County

Below is an interactive map that helps us visualize the change in population in Rutherford, Davidson, and Williamson County. The blue color represents the current estimate of the population, while the gray represents the change in population. The darker the shade, the more change or higher population estimate. Move the slider to view both.

Population Estimate & Population Change Map

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("60e19f7ff69b7368570b2a6560506673a1e4017c")

# 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