CODEWITHSUNDEEP
rusmap
M3Lba
M3Lba

Reputation: 165

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.

Upvotes: 1

Views: 161

Answers (1)

Allan Cameron
Allan Cameron

Reputation: 174586

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

Upvotes: 1

Related Questions