9

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

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.

4 Answers 4

15

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

states %>%
  group_by(group_var) %>% 
  summarize(geometry = st_union(geometry))
1
  • 1
    This is great; thank you. Since I'm more familiar with dplyr, this is much more intuitive to me.
    – Michael A
    Mar 24, 2019 at 3:49
4

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)
2
  • This is brilliant; thank you! Out of curiosity, what is the mean function actually taking the mean over if there is no non-geometric data in the data frame?
    – Michael A
    Mar 21, 2019 at 17:41
  • Running x suggests that factors get turned into NA's. Mar 21, 2019 at 18:10
3

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

0
2

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)

1
  • small typo: plot(Columbus) should be plot(columbus)!?
    – Matifou
    Jan 24 at 17:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.