Can't get regional geometries to work with maps of the US (R)

Question

TLDR: Trying to graph some simple US data by census regions and I wanted to put an outline surrounding the regions, but I can't seem to get the regional shapefiles to work.

Hi All,

I'm new to graphing geographic data in R and am trying something simple (US data by census region) but can't quite get some of my design elements to work. I was trying the code from this post (Mapping by US Census Divisions in ggplot2) and can't seem to get the shapefile for the outlines of the census regions to work properly. Here's the code from the other post:

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

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)

However, whenever I get to the block to create division polygons I get this error:

Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 236 crosses edge 239

I imagine this is something simple that I'm doing wrong with the shapefiles, but I'm not sure what the issue is, or how to solve it. I tried turning off spherical geometry and using st_make_valid() but neither of them solved the issue. How should I try to fix this?

Answer 1

Since you are already using {tigris} for getting the census divisions you can consider it also as the source for the state polygons; the package was made for that.

By doing that you would sidestep the by now superseded {maps} package. It has had a long and fruitful life - did you know it comes from S days, and is a few months older than R itself? But it does belong to a certain (less demanding) age and we really should thank it for its service, and move on.

So consider a small update of your code; what it does is:

• keeps the geometry of tigris::states() present / ie. removes the sf::st_drop_geometry() call
• filters the states by abbreviation to get rid of far flung territories such as the Virgin Islands, Guam & Alaska and Hawaii, leaving the lower 48
• the rest of your code - creating division polygons and plotting them - is left basically unchanged, only replacing the states object by the div_dat from tigris

I believe the map serves the intended purpose, and in addition the code is somewhat more concise.

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

div_dat <- states(cb = FALSE, resolution = '20m') %>%
  filter(!STUSPS %in% c('HI', 'AK', 'PR', 'GU', 'VI', 'AS', 'MP')) %>% 
  select(NAME, DIVISION) 

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

# plot it
ggplot() + 
  theme_void() +
  geom_sf(data = div_dat, 
          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

Two ways around this, and I'm not sure which is better. #1 is to change the crs (I used web mercator '3857'), make it valid, and then change the crs back to the original. #2 is to set sf_use_s2(FALSE).

Both examples below:

library(tidyverse)
library(tigris)
library(sf)
library(maps)
library(patchwork)

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)
#> Joining with `by = join_by(ID)`

# option 1change crs, make valid, change back
orig_crs <- st_crs(states)
states_change_crs <- states %>%
  st_transform(3857) %>%
  st_make_valid() %>%
  st_transform(crs = orig_crs)

# option 2 sf_use_s2(FALSE)
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
states_s2false <- states


# create division polygons
div_1 <- states_change_crs %>%
  group_by(DIVISION) %>% 
  summarize()

div_2 <- states_s2false %>%
  group_by(DIVISION) %>%
  summarize()

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

# plot #2
p2 <- ggplot() +
  theme_void() + 
  geom_sf(data = states_s2false,
          aes(fill = as.numeric(DIVISION)),
          color = 'white') +
  geom_sf(data = div_2,
          color = 'black',
          fill = NA,
          size = 1) + 
  scale_fill_viridis_c() + 
  coord_sf(crs = 5070) +
  labs(fill = NULL)

p1 + p2

^{Created on 2023-04-19 by the reprex package (v2.0.1)}