0

When using rasterToPolygons within the raster package each cell that meets the formula criteria becomes its own polygon:

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

I however want each polygon that has an adjacent side or corner to be part of one larger polygon, decreasing the number of total polygons. Is there any way to accomplish this?

Val
  • 6,585
  • 5
  • 22
  • 52
Forrest Williams
  • 125
  • 8

2 Answers2

0

OLD ANSWER:

You can use the argument dissolve=TRUE

library(raster)
r <- raster(nrow=18, ncol=36)
r[] <- sample(2, ncell(r), replace=TRUE)
pol <- rasterToPolygons(r, dissolve=TRUE)
plot(pol)

NEW ANSWER

If you do not care about the values, you can do something like this

Your example data

library(raster)
r <- raster(nrow=18, ncol=36)
r[] <- runif(ncell(r)) * 10
r[r>8] <- NA

Set all values cells you want to one value, all others to NA

x <- reclassify(r, rbind(c(-Inf, 6, NA), c(6, Inf, 1)))

pol <- rasterToPolygons(x, dissolve=TRUE)

Note that pol now has only 1 (multi-)polygon. If you want to separate the non-connected parts, you can do 

pols <- disaggregate(pol)
pols
#class       : SpatialPolygonsDataFrame 
#features    : 80

Note that diagonally adjecent polygons are separate from each other as they cannot for a valid single polygon (it would be self-intersecting).

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • The dissolve option only works here because you've forced every cell in the raster to have the same value. What I'm looking for is a way to group together all neighboring output polygons regardless of their inherited raster value. – Forrest Williams Apr 12 '18 at 18:41
  • I updated the answer now that I think that I understand your question better – Robert Hijmans Apr 12 '18 at 19:16
-1

This can be accomplished by using the poly2nb function in the spdep package to define the neighbors of each polygon, using the function created below to create a vector of region assignments, using spCbind from the maptools package to bind regions to pol, then finally dissolving 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)
Forrest Williams
  • 125
  • 8