Animated cartograms with NHGIS data

nhgis_cartogram.R

library(purrr)
library(ipumsr) # requires development version: devtools::install_github("mnpopcenter/ipumsr")
library(rmapshaper)
library(ggplot2)
library(gganimate) # devtools::install_github("thomasp85/gganimate")
library(stringr)
library(dplyr)
library(tidyr)
library(cartogram)
library(sf)

# Data and shape files from NHGIS:
# Go to nhgis.org -> Get Data
# Select the Total Population Time Series data that goes back to 1790
# Select the 2010 state shape files
# Submit extract (Make an account / log in) -> Default file options are okay
extract_number <- "0039" # Change this to match your extract

exclude_states <- c(
  "G020", # Alaska
  "G150", # Hawaii
  "G720" # Puerto Rico
)

data <- read_nhgis(paste0("data/nhgis", extract_number, "_csv.zip"), var_attrs = NULL) %>%
  filter(!GISJOIN %in% exclude_states) %>%
  select(GISJOIN, starts_with("A00AA")) %>%
  gather("year", "pop", -GISJOIN) %>%
  extract(year, "year", "A00AA([0-9]{4})$", convert = TRUE) %>%
  filter(year >= 1950 & !is.na(pop))

shapes <- read_ipums_sf(
  paste0("data/nhgis", extract_number, "_shape.zip"), 
  matches("us_state_2010"),
  vars = GISJOIN,
  verbose = FALSE
) %>%
  filter(GISJOIN %in% data$GISJOIN) %>% 
  ms_simplify(0.001)

cartograms <- inner_join(shapes, data, by = "GISJOIN") %>%
  group_by(year) %>% 
  nest(.key = "sf") %>%
  mutate(sf = map(sf, ~cartogram_cont(., "pop"))) %>%
  mutate(
    total_pop = map_dbl(sf, ~sum(.$pop)),
    total_area = map_dbl(sf, ~sum(st_area(.)))
  ) %>%
  mutate(
    cur_scalar = total_area / total_area[year == 2010],
    desired_scalar = total_pop / total_pop[year == 2010],
  ) %>%
  mutate(scalar_adj = desired_scalar / cur_scalar) %>%
  mutate(
    sf = map2(sf, scalar_adj, function(.x, .y) {
      .x$geometry <- .x$geometry * sqrt(.y)
      .x
    })
  ) %>%
  unnest() %>%
  st_as_sf() %>%
  mutate(state = paste0("Year: ", year))


area_map <- shapes %>%
  mutate(area = as.numeric(st_area(.))) %>%
  cartogram_cont("area") %>%
  mutate(
    geometry = geometry * sqrt(cartograms$scalar_adj[cartograms$year == 1950][1]),
    state = "Map"
  )

plot_data <- bind_rows(area_map, cartograms) %>%
  st_as_sf() %>% 
  mutate(state = factor(state, levels = c("Map", unique(cartograms$state))))

make_caption <- function(x) {
  if (x == "Map") return("")
  paste0("Total populaiton of Continental US: ", scales::comma(plot_data$total_pop[plot_data$state == x][1]))
}

ggplot(plot_data) + 
  geom_sf() +
  coord_sf(datum = NA) + 
  theme_minimal() + 
  transition_states(state, 2, 1) + 
  labs(
    title = "Animated population cartograms (Continental US) via NHGIS data",
    subtitle = "{closest_state}",
    caption = "{make_caption(closest_state)}"
  )