The interactive map below shows the changes in population in the greater Nashville area between 2018 and 2023. The map has a moveable filter that allows viewers to see population estimates for both years.

Code:

options(tigris_use_cache = TRUE)

# Transmitting your Census API key 

census_api_key("35989942993f23ba972685ee05578b3cdb4636bb")

# 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