#States need to be in equal area projection
#5070 is EPSG code for the USGS Alber's projection
states <- st_transform(states, crs = 5070)
#Create centroids of states
cntr <- st_centroid(states)
#Bind these coordinates to the centroid feature
cntr <- cbind(cntr, st_coordinates(cntr))
#Get centroid for conterminous US
xy <- states %>% st_union() %>%
  st_centroid() %>% st_coordinates()
#Define column 'region' and break US up into regions
cntr$region <- 'NE'
cntr$region[cntr$X > xy[1] & cntr$Y < xy[2]] <- 'SE'
cntr$region[cntr$X < xy[1] & cntr$Y < xy[2]] <- 'SW'
cntr$region[cntr$X < xy[1] & cntr$Y > xy[2]] <- 'NW'
#Region IDs back to state polygons
states$region <- cntr$region
#Dissolve by region ID
regions <- states %>%
  group_by(region) %>%
  summarise()
plot(regions$geometry)
plot(cntr['region'], pch = 19, add=T)