2

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.

4
  • plot(coordinates(out[[1]][[200]]), type="l")?
    – rcs
    Dec 1, 2013 at 22:09
  • okay, can you make a reproducible example?
    – rcs
    Dec 1, 2013 at 23:04
  • look, I made some changes in my data as well as question @rcs
    – ToNoY
    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.
    Jan 21, 2014 at 16:24

3 Answers 3

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.

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)
0

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)
0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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