Combining districts in shapefile to produce custom regions in R

Question

I am working with a dataset that uses custom regions ("North", "South" and "West") and I would need to plot them onto a map with tmap. As a first step, I assigned the different districts to the respective regions and merged them into the shapefile and then merge the data with the shapefile.

customregion <- data.frame(district = c("Area1", "Area2", 
"Area3", "Area4")

region = c("North", "North","South","South"))

shape <- sf::read_sf("./gadm40_ABC_1.shp") # read shapefile

shape <- merge(shape, customregion, by = "NAME_1")

mapdata <- merge(shape, data, by.x = "NAME_1", by.y = "district")

tm_shape(data) +
  tm_fill(midpoint = NULL, "indicator", palette="Reds", style="cont") +
  tm_text(text = "maplabel")

Not unsurprisingly, I get a map like the one below, where the regional data is plotted into each district border.

It seems I would need to aggregate the different borders in the shapefile into proper regions.

How do I do this in R?

[![enter image description here][1]][1]

  [1]: https://i.stack.imgur.com/WryPs.png

Answer 1

The geoprocess you are searching is called Dissolve. According to GIS WIKI:

Dissolve is an application of the conceptual operators that aggregates features often referred to as 'Merge' or 'Amalgamation.' It is a process in which a new map feature is created by merging adjacent polygons, lines, or regions that have a common value for a specified attribute.

This fits in your description. With sf objects, you can use dplyr functions in order to apply this geoprocess.

Some reproducible data:

library(sf)
library(dplyr)
library(rnaturalearth)

v <- ne_countries()

v <- st_as_sf(v)

v %>% filter(region_un == 'Africa') -> v

plot(v[,'subregion'])

For dissolve geometries, the first step is grouping by one (or more) attribute, then summarise geometries:

v %>% group_by(subregion) %>% summarise() %>% plot()

There are some issues with this example’s geometries. You can use sfheaders library for improving the output:

library(sfheaders)

v %>% group_by(subregion) %>% summarise() %>%
  sf_remove_holes() %>% plot()