R focal function: Calculate square of difference between raster cell and neighborhood and find the m-CodePudding

I am trying to calculate the square of the difference between a raster cell i and its neighbors js (i.e.,(j-i)^2) in a 3 x 3 neighborhood, and then calculate the mean value of those differences and assign that result to cell i. I found this answer, given by Forrest R. Stevens, that comes close to what I want to achieve, but I have only one raster (not a stack) with 136710 cells (1 089 130 combinations with the adjacent function), so a for loop is taking forever.

I want to use the function focal from the raster package, so the for loop is only run for the 3x3 matrix, but it is not working for me. Here is an example using Forrest R. Stevens' code I mentioned above:

r <- raster(matrix(1:25,nrow=5))
r[] <-c(2,3,2,3,2,
        3,2,3,2,NA,
        NA,3,2,3,2,
        NA,2,3,2,3,
        2,3,2,3,NA)
##  Calculate adjacent raster cells for each focal cell:
a <- raster::adjacent(r, cell=1:ncell(r), directions=8, sorted=T)
# Function
sq_dff<- function(w){
  ##  Create column to store calculation:
  out <- data.frame(a)
  out$sqrd_diff <- NA
  ##  Loop over all focal cells and their adjacencies,
  ##    extract the values across all layers and calculate
  ##    the squared difference, storing it in the appropriate row of
  ##    our output data.frame:
  cores <- 8
  beginCluster(cores, type='SOCK')
  for (i in 1:nrow(a)) {
    print(i)
    out$sqrd_diff[i] <- (r[a[i,2]]- r[a[i,1]])^2
    print(Sys.time())
  }
  endCluster()
  ##  Take the mean of the squared differences by focal cell ID:
  r_out_vals <- aggregate(out$sqrd_diff, by=list(out$from), FUN=mean,na.rm=T)
  names(r_out_vals)<- c('cell_numb','value')
  return(r_out_vals$value)
}


r1 <- focal(x=r, w=matrix(1,3,3), fun=sq_dff)

The function works well if I apply it like this: r1 <-sq_dff(r), and using #r_out <- r[[1]]; #r_out[] <- r_out_vals$value; return(r_out) (as suggested by. Forrest R. Stevens in his answer) instead of return(r_out_vals$value)

But, when I apply it inside the focal function as written above, it returns a raster with values for only the nine cells in the center and all of them with the same value of 0.67 assigned.

Thanks!

CodePudding user response：

You could try this:

library(terra)
r <- rast(matrix(1:25,nrow=5))
r[] <-c(2,3,2,3,2,
        3,2,3,2,NA,
        NA,3,2,3,2,
        NA,2,3,2,3,
        2,3,2,3,NA)
        
f <- function(x) {
    mean((x[-5] - x[5])^2, na.rm=TRUE)
}

rr <- focal(r, 3 ,f)
plot(rr)
text(rr, dig=2)