Joining nearest neighbor small polygons using R

Question

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

require(sp)
require(raster)
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
plot(pol100[[1]][[500]])
#if want to know the area of a particular polygon (say polygon#100)
#(all polygons are in fact of same size)
pol100@polygons[[100]]@area

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.

Answer 1

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.

library(raster)
library(spdep)
library(maptools)

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)
plot(pol)

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
    }
  }
  group
}

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

Answer 2

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({
        sf::st_as_sf(stars::st_as_stars(r),
            as_points = FALSE, merge = TRUE
        ) %>% sf::st_make_valid() # %>% sf::as_Spatial()
    })
    sf::st_crs(sf_poly) = crs
    sf_poly
}
pol100 <- rast2poly(r)

Answer 3

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

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

Solution

p <- as.polygons(x)

See

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