How do I combine geometries in a shapefile based on a grouping variable?

Question

Here is a simple example. I have a shapefile of three US states:

This is the code that creates the shapefile.

library(USAboundaries)
library(sf)
library(dplyr)

states <- us_states(map_date = "2000-01-01", resolution = "high", states = c("CA", "OR", "WA")) %>%
  mutate(group_var = if_else(state_abbr == "CA", 1, 2))
plot(st_geometry(states))

Now, I'm trying to combine the geometries based on the value of group_var, so California should stand alone while Oregon and Washington get lumped together into a single geometry. Unfortunately, st_combine doesn't take a grouping variable, and although the aggregate function in sf looks promising, code like this throws an error that group_var cannot be found.

x <- aggregate(states, group_var, mean)

Furthermore, aggregate requires a function as a third argument, presumably because aggregate in the stats package applies the function to the data, but in this case, there isn't any data to apply a function to. I'm just trying to combine/aggregate the shapefiles.

Answer 1

You could also do this using dplyr's group_by() and summarize() functions:

states %>%
  group_by(group_var) %>% 
  summarize(geometry = st_union(geometry))

Answer 2

You were close; just cast the second argument in aggregate() as a list, like this:

x <- aggregate(states, 
               by = list(states$group_var),
               FUN = mean)

x
plot(x)

Answer 3

You could also simplify Chris Merkord's answer with: states %>% group_by(group_var) %>% summarize()

Answer 4

You are basically describing a dissolve operation. This can easily be done using the rgeos package. Note that I switched datasets because I was having difficulty installing the USAboundaries and USAboundariesData packages.

library(rgeos)
library(sp)
library(spData)
library(rgdal)

columbus <- rgdal::readOGR(system.file("shapes/columbus.shp", 
                           package="spData")[1])
   print(columbus$NSB)
     plot(columbus)

diss <- rgeos::gUnaryUnion(columbus, id = columbus$NSB)
  plot(diss, col=c("red","green"))
    plot(columbus, add=TRUE)

For a more modern interface using sf the equivalent is sf::st_union. Keep in mind that the default precision can have an effect on how internal boundaries are resolved. This can be changed using something like: sf::st_set_precision(1e8)