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.

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

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.