How start with one municipality. Find the first bordering ring - defined as the set of municipalities bordering the first. The add the second ring, all those municipalities bordering a municipality in the first ring. Then …. third ring … then ….

wd <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//R_space_examples//municipality_rings//"
setwd(wd)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
mun <- read_sf("mun_shape.shp")
mun <- mun[order(mun$kode),]
out <- st_intersects(mun)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
plot(mun[out[[1]],])

# Find the first neighbours 
out <- st_intersects(mun)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
i_1 <-  unique(unlist(out[[1]]))
m_1 <- mun$kode[i_1]
plot(mun[i_1,])

# Find the second neighbours ( union with first neighbours )
i_2 <- unique(unlist(out[i_1]))
m_2 <- mun$kode[i_2]
plot(mun[i_2[which(!i_2%in%i_1)],])

plot(mun[i_2,])

# Find the third set
i_3 <- unique(unlist(out[i_2]))
m_3 <- mun$kode[i_3]
plot(mun[i_3,])

# Find the third set
i_4 <- unique(unlist(out[i_3]))
m_4 <- mun$kode[i_4]
plot(mun[i_4,])