Merging small polygons with largest neighbour in R

Question

I have a bunch of polygons representing land cover classes, however I don't need this much detail. I would like to merge small polygons (ie < 150m2) with it's largest neighbour. This would be similar to "Eliminate" in ArcGIS Pro or "eliminate selected polygons" in QGIS. This process will be repeated so I'd like to script it in R and save some time.

The only way I could think to do it was to add a column indicating which polygons were under 150m2, then somehow merging them to nearest neighbor.

#Packages
library(terra)
library(sf)
library(dplyr)

#Adding polygons from Terra and converting to sf
v <- vect(system.file("ex/lux.shp", package="terra"))
v <- st_as_sf(v[c(1:12)])


#Adding a column indicating which polygons are under 150m2
mutate(v, area_thresh = case_when(AREA < 150 ~ 1,
                               AREA > 150 ~ 0))

Answer 1

This is an interesting problem, as you need to iterate over rows of changing shapefile (your number of rows will be decreasing as you merge smaller objects).

For a possible approach consider this piece of code that builds upon the well known & much loved NC shapefile that ships with {sf}.

Note that I am using the county_id as a unique identifier; a merged polygon will keep ID of the larger county.

library(sf)
library(dplyr)

# the one & only NC shapefile... as geometry only
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  select(CNTY_ID)

# limit for merging - in this case median area
limit <- median(st_area(shape))

# initialize iterator
i <- 1

# iterate over rows of shapefile
while (i < nrow(shape)) {
  
  # check if area over or under limit
  if(st_area(shape[i, ]) < limit) {
    
    # find id of largest neighbor
    largest_neighbor <- shape[unlist(st_touches(shape[i, ], shape)), ] %>% 
      mutate(area = st_area(.)) %>% 
      top_n(1, area) %>% 
      pull(CNTY_ID)
    
    # merge offending shape with largest neighbor
    wrk_shape <- st_union(st_geometry(shape[i, ]), 
                          st_geometry(shape[shape$CNTY_ID == largest_neighbor, ])) %>% 
      st_make_valid() %>% 
      st_as_sf() %>% 
      mutate(CNTY_ID = largest_neighbor)
    
    # rename geometry, column to enable binding
    st_geometry(wrk_shape) = attr(shape, "sf_column")
    
    # remove offending shape and unmerged neighbor
    shape <- shape[-i, ]
    shape <- filter(shape, !CNTY_ID == largest_neighbor) 
    
    # add merged shape back in
    shape <- bind_rows(shape, wrk_shape)
    
  }
  
  # don't forget to update iterator :)
  i <- i + 1
    
}

plot(st_geometry(shape))

Answer 2

You can try terra::combineGeoms.

Example data

library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))

Solution

th <- v$AREA < 150
x <- combineGeoms(v[!th], v[th])

Illustration

v$th <- th
plot(v, "th")
lines(x, lwd=3, col="blue")