[R-sig-Geo] groupping contiguous identical raster cells

[Babak Naimi]

Dear Robert and list,

I wrote a function, named 'group', to find the spatially contiguous groups of identical valued pixels in a raster layer and assign them unique integer identifiers (same as the GROUP function in Idrisi). I used the Rasterlayer object in the function to get and set values in the grouping procedure. Unfortunately the procedure is not efficient. When I converted the raster into matrix (in 'group2' function), the efficiency is incredibly improved but I think it still has problem with bigger rasters. It is not also a memory-safe function!

I was wondering if anyone can help to modify the function and enhance the efficiency. Following I show the results of some experiments on these functions (the source of the functions is attached at the end).

r1 <- c(rep(1,100),rep(2,100),sample(c(3,4),200,replace=T))
r1 <- raster(matrix(r1,nc=20,nr=20,byrow=T))

system.time(
  g1 <- group(r1)
  )

   user  system elapsed 
   2.56    0.00    2.59

system.time(
  g2 <- group2(r1)
  )

  user  system elapsed 
   0.06    0.00    0.06


# bigger raster:
r2 <- c(rep(1,10000),rep(2,10000),sample(c(3,4),20000,replace=T))
r2 <- raster(matrix(r2,nc=200,nr=200,byrow=T))

system.time(
  g3 <- group(r2)
  )

   user  system elapsed 
1488.39    3.54 1547.00

system.time(
  g4 <- group2(r2)
  )
  user  system elapsed 
 230.47    0.52  237.54

#--------------------
# The source of functions:
cellNeighbors<-function(cell,nr,nc,type="queen") {
  row <- trunc((cell-1)/nc) + 1
  col <- cell - ((row-1) * nc)
  
  if (type == "rook") {
    r <- c(-1,0,1,0)
    c <- c(0,1,0,-1)
  } else {
    r <- c(-1,-1,-1,0,1,1,1,0)
    c <- c(-1,0,1,1,1,0,-1,-1)
  }
  row <- row - r
  col <- col - c
  row[row < 1 | row > nr] <- NA
  col[col < 1 | col > nc] <- NA
  rc <- cbind(row, col)
  v <- (rc[,1]-1) * nc + rc[,2]
  w <- which(is.na(v)) 
  if (length(w) > 0) v <- v[-w]
  v
}
#---------
group <- function(x,type="queen") {
  g <- raster(x)
  nc <- x at ncols
  nr <- x at nrows
  cells <- 1:ncell(x)
  w <- which(is.na(x[cells]))
  if (length(w) > 0) cells <- cells[-w]
  groupID <- 0
  Neigh.cells <- c()
  while (length(cells) > 0) {
    groupID <- groupID + 1
    cell <- cells[1]
    g[cell] <- groupID
    cells <- cells[-1]
    cell.n <- cellNeighbors(cell,nr,nc,type=type)
    w <- which(x[cell.n] == x[cell])
    if (length(w) > 0) Neigh.cells <- c(Neigh.cells,cell.n[w])
    while (length(Neigh.cells) > 0) {
      cell <- Neigh.cells[1]
      Neigh.cells <- Neigh.cells[-1]
      w <- which(cells == cell)
      if (length(w) > 0) {
        cells <- cells[-w]
        g[cell] <- groupID
        cell.n <- cellNeighbors(cell,nr,nc,type=type)
        ww <- which(x[cell.n] == x[cell])
        if (length(ww) > 0) Neigh.cells <- c(Neigh.cells,cell.n[ww])
      }
    }
  }
  g
}

#--------
group2 <- function(x,type="queen") {
  gr <- raster(x)
  nc <- x at ncols
  nr <- x at nrows
  cells <- 1:ncell(x)
  x <- t(as.matrix(x))
  g <- x
  g[] <- NA
  w <- which(is.na(x[cells]))
  if (length(w) > 0) cells <- cells[-w]
  groupID <- 0
  Neigh.cells <- c()
  while (length(cells) > 0) {
    groupID <- groupID + 1
    cell <- cells[1]
    g[cell] <- groupID
    cells <- cells[-1]
    cell.n <- cellNeighbors(cell,nr,nc,type=type)
    w <- which(x[cell.n] == x[cell])
    if (length(w) > 0) Neigh.cells <- c(Neigh.cells,cell.n[w])
    while (length(Neigh.cells) > 0) {
      cell <- Neigh.cells[1]
      Neigh.cells <- Neigh.cells[-1]
      w <- which(cells == cell)
      if (length(w) > 0) {
        cells <- cells[-w]
        g[cell] <- groupID
        cell.n <- cellNeighbors(cell,nr,nc,type=type)
        ww <- which(x[cell.n] == x[cell])
        if (length(ww) > 0) Neigh.cells <- c(Neigh.cells,cell.n[ww])
      }
    }
  }
  gr[] <- t(g)
  gr
}
#---------------------

Best regards,
Babak

Babak Naimi
PhD Candidate, Natural Resources Department,
Faculty of Geo-Information Science and Earth Observations (ITC), University of Twente, 
P.O. Box 217, 7500 AE Enschede, The Netherlands
Tel: +31 53 4874212