Get values from adjacent cells raster R

Question

I need to extract values from cells adjacent to the my focus cells in a raster and I need to be able to link these values with the focus cell. So far I can extract the values from the adjacent cells but the outcome is not in an organised form.

 library(raster)
 r <- raster(matrix(runif(100), 10))
 cells <- c(34,22,50,10)
 cells_ad <- adjacent(r,cells, directions = 8, pairs=FALSE)
 extract (r,cells_ad)

I need to have a data frame with the values extracted for each of the adjacent cells with each column being one adjacent cell (ideally I would know the order - south, southwest, west, northwest...) and each row a focus cell. Following my example above the data frame would have 8 columns and 4 rows.

Need to do this for a massive raster file and millions of points to be extracted. So computation time should be considered here.

I appreciate any help.

Answer 1

Here is a simple Queen-case adjacent function that includes cells outside the raster (as NA).

adj <- function(x, rr, ngb, global) {
    a <- x + ngb
    # the below line is not strictly needed, so it should be faster without it
    # a[a < 1 | a > ncell(rr)] <- NA
    if (!global) {
        col <- colFromCell(rr, x)
        if (col == 1) {
            a[c(1,4,6)] <- NA
        } else if (col==nc) {
            a[c(3,5,8)] <- NA
        }
    }
    a
}

Your example data

library(raster)
r <- raster(matrix(runif(100), 10))
cells <- c(34,22,50,10)

Apply the function

nc <- ncol(r)
nbrs <- c(-nc-1, -nc, -nc+1, -1, +1, +nc-1, +nc, +nc+1) 
ad <- t(sapply(cells, adj, rr=raster(r), ngb=nbrs, global=FALSE))   
colnames(ad) <- c("NW", "N", "NE", "W", "E", "SW", "S", "SE")

Extract the values

ad[] <- extract(r, as.vector(ad))

Answer 2

See below my (sloppy) solution. I noticed the adjacent() function omit NA cells, so I modified it to remove that. The orientation for the extracted cells follows this order: "nw","w","sw","ne","e","se","n","s".

library(raster)
r <- raster(matrix(runif(100), 10))
cells <- c(34,22,50,10)

# modified adjacent function to exclud the na.omit() part, therefore it returns all adjacent cells around the cell
adjacent_mod = function (x, cells, directions = 4, pairs = TRUE, target = NULL, 
                         sorted = FALSE, include = FALSE, id = FALSE) 
{
  .isGlobalLonLat <- function(x) {
    res <- FALSE
    tolerance <- 0.1
    scale <- xres(x)
    if (isTRUE(all.equal(xmin(x), -180, tolerance=tolerance, scale=scale)) & 
        isTRUE(all.equal(xmax(x),  180, tolerance=tolerance, scale=scale))) {
      if (couldBeLonLat(x, warnings=FALSE)) {
        res <- TRUE
      }
    }
    res
  }


  if (is.character(directions)) {
    directions <- tolower(directions)
  }
  x <- raster(x)
  r <- res(x)
  xy <- xyFromCell(x, cells)
  mat <- FALSE
  if (is.matrix(directions)) {
    stopifnot(length(which(directions == 0)) == 1)
    stopifnot(length(which(directions == 1)) > 0)
    d <- .adjacentUD(x, cells, directions, include)
    directions <- sum(directions == 1, na.rm = TRUE)
    mat <- TRUE
  }
  else if (directions == 4) {
    if (include) {
      d <- c(xy[, 1], xy[, 1] - r[1], xy[, 1] + r[1], 
             xy[, 1], xy[, 1], xy[, 2], xy[, 2], xy[, 2], 
             xy[, 2] + r[2], xy[, 2] - r[2])
    }
    else {
      d <- c(xy[, 1] - r[1], xy[, 1] + r[1], xy[, 1], 
             xy[, 1], xy[, 2], xy[, 2], xy[, 2] + r[2], xy[, 
                                                           2] - r[2])
    }
  }
  else if (directions == 8) {
    if (include) {
      d <- c(xy[, 1], rep(xy[, 1] - r[1], 3), rep(xy[, 
                                                     1] + r[1], 3), xy[, 1], xy[, 1], xy[, 2], rep(c(xy[, 
                                                                                                        2] + r[2], xy[, 2], xy[, 2] - r[2]), 2), xy[, 
                                                                                                                                                    2] + r[2], xy[, 2] - r[2])
    }
    else {
      d <- c(rep(xy[, 1] - r[1], 3), rep(xy[, 1] + r[1], 
                                         3), xy[, 1], xy[, 1], rep(c(xy[, 2] + r[2], 
                                                                     xy[, 2], xy[, 2] - r[2]), 2), xy[, 2] + r[2], 
             xy[, 2] - r[2])
    }
  }
  else if (directions == 16) {
    r2 <- r * 2
    if (include) {
      d <- c(xy[, 1], rep(xy[, 1] - r2[1], 2), rep(xy[, 
                                                      1] + r2[1], 2), rep(xy[, 1] - r[1], 5), rep(xy[, 
                                                                                                     1] + r[1], 5), xy[, 1], xy[, 1], xy[, 2], rep(c(xy[, 
                                                                                                                                                        2] + r[2], xy[, 2] - r[2]), 2), rep(c(xy[, 2] + 
                                                                                                                                                                                                r2[2], xy[, 2] + r[2], xy[, 2], xy[, 2] - r[2], 
                                                                                                                                                                                              xy[, 2] - r2[2]), 2), xy[, 2] + r[2], xy[, 2] - 
               r[2])
    }
    else {
      d <- c(rep(xy[, 1] - r2[1], 2), rep(xy[, 1] + r2[1], 
                                          2), rep(xy[, 1] - r[1], 5), rep(xy[, 1] + r[1], 
                                                                          5), xy[, 1], xy[, 1], rep(c(xy[, 2] + r[2], 
                                                                                                      xy[, 2] - r[2]), 2), rep(c(xy[, 2] + r2[2], 
                                                                                                                                 xy[, 2] + r[2], xy[, 2], xy[, 2] - r[2], xy[, 
                                                                                                                                                                             2] - r2[2]), 2), xy[, 2] + r[2], xy[, 2] - 
               r[2])
    }
  }
  else if (directions == "bishop") {
    if (include) {
      d <- c(xy[, 1], rep(xy[, 1] - r[1], 2), rep(xy[, 
                                                     1] + r[1], 2), xy[, 2], rep(c(xy[, 2] + r[2], 
                                                                                   xy[, 2] - r[2]), 2))
    }
    else {
      d <- c(rep(xy[, 1] - r[1], 2), rep(xy[, 1] + r[1], 
                                         2), rep(c(xy[, 2] + r[2], xy[, 2] - r[2]), 2))
    }
    directions <- 4
  }
  else {
    stop("directions should be one of: 4, 8, 16, \"bishop\", or a matrix")
  }
  if (include) 
    directions <- directions + 1
  d <- matrix(d, ncol = 2)
  if (.isGlobalLonLat(x)) {
    d[, 1] <- (d[, 1] + 180)%%360 - 180
  }
  if (pairs) {
    if (mat) {
      cell <- rep(cells, each = directions)
    }
    else {
      cell <- rep(cells, directions)
    }
    if (id) {
      if (mat) {
        ID <- rep(1:length(cells), each = directions)
      }
      else {
        ID <- rep(1:length(cells), directions)
      }
      d <- cbind(ID, cell, cellFromXY(x, 
                                      d))
      attr(d, "na.action") <- NULL
      colnames(d) <- c("id", "from", "to")
      if (!is.null(target)) {
        d <- d[d[, 3] %in% target, ]
      }
    }
    else {
      d <- cbind(cell, cellFromXY(x, d))
      attr(d, "na.action") <- NULL
      colnames(d) <- c("from", "to")
      if (!is.null(target)) {
        d <- d[d[, 2] %in% target, ]
      }
    }
    if (sorted) {
      d <- d[order(d[, 1], d[, 2]), ]
    }
  }
  else {
    d <- as.vector(unique(cellFromXY(x, d)))
    if (!is.null(target)) {
      d <- intersect(d, target)
    }
    if (sorted) {
      d <- sort(d)
    }
  }
  d
}

# extract adjacent cell numbers
cells_ad <- adjacent_mod(r,cells, directions = 8, pairs=T, sorted=T, id=F)

# table them into a data.frame with the orientation
i=1
for (i in 1:length(unique(cells))) {
  if (i == 1) {
    df = data.frame(matrix(c(cells[i], cells_ad[which(cells_ad[,1]==cells[i]),2]), ncol=9))
  } else {
    df = rbind(df, data.frame(matrix(c(cells[i], cells_ad[which(cells_ad[,1]==cells[i]),2]), ncol=9)))
  }
}
names(df) = c("cell","nw","w","sw","ne","e","se","n","s")

# extract values from raster r
i=1
df_values = df
for (i in 1:dim(df)[1]) {
  df_values[i,] = extract (r, as.numeric(df[i,]))
}