Search code examples
rusmap

Error message when converting usmap object to sf object


I am attempting to create a map with county level data in several US states. Since I'm including Alaska, I cannot use the maps::map() function and have elected to use the usmap::us_map() function. The result is a df that produces a lovely map.

df<-usmap::us_map('counties')%>%
  filter(abbr %in% c('AK','CO','ID','MT','UT','OR','WA','WY','ND','SD'))%>%
  rename(state=abbr)
usmap::plot_usmap(regions='counties', include=c('AK','CO','ID','MT','UT','OR','WA','WY','ND','SD'))

enter image description here

However, the plotting options that exist within the usmap package aren't quite what I need. I want to convert it to an sf object. I used the code below,

MtWc <- lapply(split(df, list(df$county,df$piece),drop=T), function(x) st_polygon(list(cbind(x$x, x$y))))
MtW  <- st_sfc(MtWc, crs = usmap_crs()@projargs)
MtW  <- st_sf(data.frame(fips = unique(d$fips), county = names(MtWc), geometry = MtW))

and although it works for some states(eg Colorado only), it does not work for all and results in the following error:

Error in MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) : polygons not (all) closed

I've tried a variety of methods to avoid this, including

df<-cleangeo::clgeo_Clean(df)

Basically, does anyone know a) how to get an sf object with county-level polygons for the listed region or b) how to close the polygons to make this work?

FYI, I chose to split by county and piece due to the states that have multiple pieces of each county, usually islands.


Solution

  • You can split the data frame by state, county and piece, then create an sf for each of these and bind them together:

    library(sf)
    
    df <- usmap::us_map('counties') %>%
      rename(state = abbr) %>%
      filter(state %in% c('AK','CO','ID','MT','UT','OR','WA','WY','ND','SD'))
    
    dat <- do.call('rbind', split(df, paste(df$state, df$county, df$piece)) |>
      lapply(function(x) st_sf(county = x$county[1], state = x$state[1],
                               st_sfc(st_polygon(list(cbind(x$x, x$y)))), 
                               crs = usmap::usmap_crs())))
    

    Now dat is an sf data frame:

    library(ggplot2)
    
    ggplot(dat) + geom_sf(fill = 'white')
    

    enter image description here