Mapping by US Census Divisions in ggplot2

Question

I would like to map data by US census regions. In short, I have state-level data, which I joined to latitude/longitude coordinates using the maps package. Here are the first 20 lines of the resulting dataframe:

 df_F4 %>% select(state, division, percap_rheum, group, long, lat) %>% print(n=20)
# A tibble: 15,539 x 6
# Groups:   division [9]
   state   division           percap_rheum group  long   lat
   <chr>   <chr>                     <dbl> <dbl> <dbl> <dbl>
 1 alabama east south central       70255.     1 -87.5  30.4
 2 alabama east south central       70255.     1 -87.5  30.4
 3 alabama east south central       70255.     1 -87.5  30.4
 4 alabama east south central       70255.     1 -87.5  30.3
 5 alabama east south central       70255.     1 -87.6  30.3
 6 alabama east south central       70255.     1 -87.6  30.3
 7 alabama east south central       70255.     1 -87.6  30.3
 8 alabama east south central       70255.     1 -87.6  30.3
 9 alabama east south central       70255.     1 -87.7  30.3
10 alabama east south central       70255.     1 -87.8  30.3
11 alabama east south central       70255.     1 -87.9  30.2
12 alabama east south central       70255.     1 -87.9  30.2
13 alabama east south central       70255.     1 -88.0  30.2
14 alabama east south central       70255.     1 -88.0  30.2
15 alabama east south central       70255.     1 -88.0  30.3
16 alabama east south central       70255.     1 -88.0  30.3
17 alabama east south central       70255.     1 -88.0  30.3
18 alabama east south central       70255.     1 -88.0  30.3
19 alabama east south central       70255.     1 -87.9  30.3
20 alabama east south central       70255.     1 -87.8  30.3
# … with 15,519 more rows

I want to graph it with ggplot2, organized by census regions. I have aggregated the data by regions and can graph it at a state level as such:

graph_theme <-  theme_light() + 
  theme(
    text = element_text(size=10),
    panel.grid.major.x = element_blank(), 
    panel.grid.minor = element_blank(), 
    plot.margin = unit(c(0.75, 0.25, 0.5, 0.5), "cm")) #top, R, bottom, L) 
  
map_theme <- theme(
    axis.title.x = element_blank(), 
    axis.title.y = element_blank(), 
    axis.text.x = element_blank(), 
    axis.text.y = element_blank(), 
    axis.ticks = element_blank(), 
    panel.grid = element_blank(), 
    legend.title = element_blank(), 
    legend.position = c(0.92, 0.25), # h / v 
    legend.background = element_blank()
)

df_F4 %>%
  ggplot(aes(
    long, 
    lat, 
    group = group)) +
  geom_polygon(aes(fill = percap_rheum), color = "white") + 
  scale_fill_viridis_c(labels = dollar_format(big.mark = ","), direction = -1) + 
  graph_theme + 
  map_theme

Which results in this picture:

You can tell that there are census divisions there, but I want to highlight them or outline them somehow. Any advice would be greatly appreciated!

Answer 1

One approach would be to create a polygon for each of the divisions, then overlay those over the state data.

An example using the data from tigris::state(). The steps are:

1. download the data and filter to CONUS,
2. create polygons of divisions,
3. plot state data with boundaries of divisions overlaid.

I also changed to a geographic crs which gives the US a bit of a curved look. Don't have to do that though.

library(tidyverse)
library(tigris)
library(sf)

# Download state data and filter states
sts <- states() %>%
  filter(!STUSPS %in% c('HI', 'AK', 'PR', 'GU', 'VI', 'AS', 'MP'))

# Summarize to DIVISION polygons, see sf::st_union
div <- sts %>%
  group_by(DIVISION) %>% 
  summarize()

# Plot it
ggplot() + 
  theme_void() +
  geom_sf(data = sts, 
          aes(fill = as.numeric(DIVISION)), 
          color = 'white') +
  geom_sf(data = div, 
          color = 'black', 
          fill = NA,
          size = 1) +
  scale_fill_viridis_c() +
  coord_sf(crs = 5070) +
  labs(fill = NULL)
  

---

EDIT: updated to use maps package. Found a useful hint here

# get data specifying which states are in which division
div_dat <- states(cb = FALSE, resolution = '20m') %>%
  st_drop_geometry() %>%
  select(NAME, DIVISION) %>%
  mutate(ID = tolower(NAME))

# get state data, convert to sf, join with division data
states <- maps::map("state", plot = FALSE, fill = TRUE) %>%
  st_as_sf() %>%
  left_join(div_dat)

# create division polygons
div <- states %>%
  group_by(DIVISION) %>% 
  summarize()

# plot it
ggplot() + 
  theme_void() +
  geom_sf(data = states, 
          aes(fill = as.numeric(DIVISION)), 
          color = 'white') +
  geom_sf(data = div, 
          color = 'black', 
          fill = NA,
          size = 1) +
  scale_fill_viridis_c() +
  coord_sf(crs = 5070) +
  labs(fill = NULL)

Answer 2

Cred to @nniloc above for the answer to this. If anyone else is following along, we left it off with some resolution problems using the geom_sf. I tried for awhile to get geom_polygon and the "maps" package data to work. It's just a higher resolution because there are more data points I think.

Anyway, I finally settled on a compromise. Can't bend the map but the state lines themselves look good. In short, I used the division overlay described above on top of the map I had previously created. It winds up looking good enough, figured I would share if it helps others.

# get data specifying which states are in which division
div_dat <- states(cb = FALSE, resolution = '20m') %>%
  st_drop_geometry() %>%
  select(NAME, DIVISION) %>%
  mutate(ID = tolower(NAME))

# get state data, convert to sf, join with division data
states <- maps::map("state", plot = FALSE, fill = TRUE) %>%
  st_as_sf() %>%
  left_join(div_dat)

# create division polygons
div <- states %>%
  group_by(DIVISION) %>% 
  summarize()

# Plot percapita spending 

ggplot() + 
  graph_theme + 
  map_theme + 
  geom_polygon(data = df_F4, 
               aes(long, lat, group = group, fill = percap_rheum), 
               color = "white") +
  geom_sf(data = div, 
          color = "#838383", 
          fill = NA,
          size = 1) + 
  scale_fill_viridis_c(labels = dollar_format(big.mark = ","), direction = -1)