The figure shows a plot of numerous polygons (zoom in) made from a raster file by

pol100 <- rasterToPolygons(r, fun=function(x){x==1}, dissolve=FALSE)

There polygon has 2059 features

> pol100
class       : SpatialPolygonsDataFrame 
nfeatures   : 2059 
extent      : 254358.6, 268808.1, 2619318, 2642087  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=46 +datum=WGS84 +units=m +no_defs 
nvariables  : 1
names       : value 
min values  :   1 
max values  :   1 

Each of the small polygons can be explored/plotted individually

#if want to plot a particular polygon
#if want to know the area of a particular polygon (say polygon#100)
#(all polygons are in fact of same size)

I want to know if there is any automated way to merge the nearest neighbor polygons? Looking at the figure tells that the process will create approximately 15-20 large polygons, for instance. Thanks in advance.

  • plot(coordinates(out[[1]][[200]]), type="l")?
    – rcs
    Commented Dec 1, 2013 at 22:09
  • okay, can you make a reproducible example?
    – rcs
    Commented Dec 1, 2013 at 23:04
  • look, I made some changes in my data as well as question @rcs
    – ToNoY
    Commented Dec 2, 2013 at 5:48
  • If it doesn't have to be the nearest neighbors, but if you only want to reduce the complexity, I think you can use, unionSpatialPolygons {maptools} and the ID= argument to base spatial union on... I did something similar in the answer of this question; gis.stackexchange.com/questions/82667/…
    – jO.
    Commented Jan 21, 2014 at 16:24

3 Answers 3


If you want to combine all polygons that neighbor eachother you can use this method. First, use the poly2nb function in the spdep package to define the neighbors of each polygon, then use the function defined below to create a vector of region assignments, next use spCbind from the maptools package to bind regions to pol, and finally dissolve over regions using the unionSpatialPolygons function from maptools. The basic structure of the created function is if at least one of the polygon's neighbors has been assigned to a group then assign polygon and neighbors to that group else assign polygon and neighbors to new group.


r <- raster(nrow=18, ncol=36)
r[] <- runif(ncell(r)) * 10
r[r>8] <- NA
pol <- rasterToPolygons(r, fun=function(x){x>6}, dissolve = T)

nb <- poly2nb(pol)

create_regions <- function(data) {
  group <- rep(NA, length(data))
  group_val <- 0
  while(NA %in% group) {
    index <- min(which(is.na(group)))
    nb <- unlist(data[index])
    nb_value <- group[nb]
    is_na <- is.na(nb_value)
    if(sum(!is_na) != 0){
      prev_group <- nb_value[!is_na][1]
      group[index] <- prev_group
      group[nb[is_na]] <- prev_group
    } else {
      group_val <- group_val + 1
      group[index] <- group_val
      group[nb] <- group_val

region <- create_regions(nb)
pol_rgn <- spCbind(pol, region)
pol2 <- unionSpatialPolygons(pol_rgn, region)

There is a more elegant way, with the help of sf and stars:

#' rast2poly
#' @param r A raster object
#' @export
rast2poly <- function(r, crs = 4326) {
    sf_poly = suppressWarnings({
            as_points = FALSE, merge = TRUE
        ) %>% sf::st_make_valid() # %>% sf::as_Spatial()
    sf::st_crs(sf_poly) = crs
pol100 <- rast2poly(r)

If you were to join the nearest neighbors of a raster cells, you would end up with just a single polygon (the extent). I understand the goal to be to combine raster cells that have the same values into single polygons. For that, you can use terra::as.polygons.

Example data

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
x <- classify(r, 3)


p <- as.polygons(x)


plot(p, "elevation")
plot(x); lines(p)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.