Drawing state boundaries in county map in recent version

Question

Using usmap, I have made a plot of values by US counties from some values stored in a variable county. Before I recently updated R and all my packages, I could easily create a county map with thick state borders and with state abbreviations using the following code:

library(ggplot2)
library(usmap)
centroid_labels = usmapdata::centroid_labels("states")
plotdf = data.frame(cbind(county$FIPS,county$SES)); plotdf[,2] = as.numeric(as.character(plotdf[,2])); names(plotdf) = c("fips", "SES");
plot_usmap(data = plotdf, values = "SES", color = NA) +  # color = NA removes county boundary lines
geom_polygon(data = usmapdata::us_map(regions = "states"), aes(x, y, group = group), fill = NA, linewidth = 1, color = "black") + # Keep state boundary lines
geom_text(data = centroid_labels, aes(x = x, y = y, label = abbr), fontface = "bold", colour = "gray12", size = 3)

With a few more lines to change colors and such, that would produce a picture like this:

However, as mentioned, since I updated packages, this doesn't work for me anymore. Specifically, it gives me and error that the object x wasn't found. Searching for an answer, it seems like others have stumbled into a similar issue. This means that I can draw the counties, but the code lines that creates state lines and state abbreviations no longer work. Removing what doesn't work, it will produce something like this:

Question: Is there any other way to draw state boundaries in this county map and also add state abbreviations?

Answer 1

Here is another approach that is closer to the original question asked.

library(ggplot2)
library(usmap)

# Extract data for states and counties
states_data <- us_map("states")
counties_data <- us_map("counties")
counties_data$SES <- rnorm(3144, 100, 100)
centroid_labels = usmapdata::centroid_labels("states")

# Map
plot_usmap("counties", color = 'NA') +
  geom_sf(data = counties_data,
          mapping = aes(geometry = geom, fill = SES),
          color = NA) +
  geom_sf(data = states_data,
          mapping = aes(geometry = geom),
          color = 'black',
          linewidth = 1,
          fill = NA) +
  geom_sf_text(data = centroid_labels,
               mapping = aes(geometry = geom, label = abbr),
               fontface = 'bold',
               color = 'gray12',
               size = 3) +
  scale_fill_gradient2(low = "blue", high = "red", mid = 'white',
                       midpoint = 0) +
  theme(legend.position = 'none')

^{Created on 2024-05-24 with reprex v2.1.0}

Answer 2

Here is one approach using the packages geodata and sf.

library(tidyverse)
library(geodata)
library(sf)

# USA states
usa_shp <- gadm(country = 'USA', level = 1, path=tempdir()) %>%
  st_as_sf()

# USA counties
usa_county <- gadm(country = 'USA', level = 2, path=tempdir()) %>%
  st_as_sf() %>% 
  filter(ENGTYPE_2 != 'Water body') %>%
  select(COUNTRY,NAME_1,NAME_2) %>%
  rename(country = COUNTRY,
         state = NAME_1,
         county = NAME_2) %>%
  mutate(SES = round(rnorm(3137, 100, 50)))

# Centroids of USA states
centroids_usa <- usa_shp %>% 
  st_centroid(of_largest_polygon = TRUE) %>%
  select(ISO_1) %>%
  separate(ISO_1,
           c('country','state_abb'))
#> Warning: st_centroid assumes attributes are constant over geometries

# Map
ggplot() +
  geom_sf(data = usa_county,
          mapping = aes(geometry = geometry, fill = SES),
          color = NA) +
  geom_sf(data = usa_shp,
          mapping = aes(geometry = geometry),
          color = "black",
          linewidth = 1,
          fill = NA) +
  geom_sf_text(data = centroids_usa,
          mapping = aes(geometry = geometry, label = state_abb),
          fontface = 'bold',
          color = 'gray12',
          size = 3) +
  coord_sf(xlim = c(-125.5, -66), ylim = c(24.3, 49.5), expand = F) +
  scale_fill_gradient2(low = "blue", high = "red", mid = 'white',
                       midpoint = 0) +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(text = element_text(size = 16),
        axis.text.x = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"),
        legend.position = 'none')
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

^{Created on 2024-05-24 with reprex v2.1.0}