I have struggled in the past to find a reliable way to merge polygons of administrative regions in the workflow of sf package. I kept encountering two problems: invalid geometries and artefacts (slivers) at the place of former boundaries.

While both issues had a rational explanation, both were rather annoying. I created a workflow that addressed the issues, but it was sort of clunky.

As I was encountering the issue of creating arbitrary regions repeatedly I decided to wrap the code in a function. It is now a part of the RCzechia package.

To demonstrate the ease of use of the function I will aggregate the 77 Czech LAU1 regions into just two multipolygons - using as key a made up column classifying the LAU1 regions into odd and even ones, based on the number of letters in the region’s name.

The nchar() function and %% modulo operator split into one bucket the odd districts - the likes of Praha - and into second one the even lettered ones, such as the likes of Brno.

library(tidyverse)
library(RCzechia)
library(sf)


poly <- okresy("low") %>% # Czech LAU1 regions as sf data frame
  mutate(oddeven = ifelse(nchar(NAZ_LAU1) %% 2 == 1, "odd", "even" )) %>% # odd or even?
  union_sf("oddeven") # ... et facta est lux

plot(poly, key.pos = 1)

The function union_sf takes as required arguments a sf data frame (in the first position, so it is pipe friendly) and a name of a column.

Two optional arguments - tolerance and planarCRS - are used in the process of buffering the original polygons to avoid slivers. The default values - tolerance of 1 meter and planar CRS EPSG:5514 - should cover majority of situations in the Czech Republic and near abroad without causing noticeable distortion.