How to create correct spatial lag variables for a raster in R?

Question

lag.listw creates incorrect spatial lag values when I use the function spdep::cell2nb. I have a raster file and want to create a new raster where each cell has the average value of its neighboring cells (spatial lag value).

The code below creates

1. a new raster from scratch and
2. calculates the neighbor matrix with cell2nb.
3. nb2listw constructs the weights that accord to each neighbor value.
4. lag.listw creates the vector of neighbor values
5. Finally I use this vector to create a new raster.

Code:

library(raster)
library(spdep)
##raster
r<-raster(nrows=7, ncols=8)

##raster values
v<-rep(0,ncell(r))
i<-sample(1:ncell(r),1)
v[i]<-1
values(r)<-v
plot(r)

##neighbor values
#neighbor list
nb<-cell2nb(nrow=nrow(r),ncol=ncol(r),type="queen")
#spatial weights matrix
nb.w<-nb2listw(nb,style="W", zero.policy=T)
#lagged values
nb.v<-lag.listw(nb.w,values(r),zero.policy=T,NAOK=T)

##new raster
nb.r<-r
values(nb.r)<-nb.v
plot(nb.r)

The first raster looks like:

The new raster with neighbor values is:

Comparing both images it becomes clear that this method values are misplaced and wrong.

The above code works only if the given raster/ cell-matrix has equal numbers of rows and columns. Test:

##raster
r<-raster(nrows=8, ncols=8)

##raster values
v<-rep(0,ncell(r))
i<-sample(1:ncell(r),1)
v[i]<-1
values(r)<-v
plot(r)

##neighbor values
#neighbor list
nb<-cell2nb(nrow=nrow(r),ncol=ncol(r),type="queen")
#spatial weights matrix
nb.w<-nb2listw(nb,style="W", zero.policy=T)
#lagged values
nb.v<-lag.listw(nb.w,values(r),zero.policy=T,NAOK=T)

##new raster
nb.r<-r
values(nb.r)<-nb.v
plot(nb.r)

Answer 1

Use the focal function from the raster package itself. It statistics based on the neighboring cell values and the own cell value. To exclude the own cell value from this calculation you have to 1) attribute a zero weight to it and 2) adapt the mean function to have one observation less.

##create base raster
r<-raster(nrows=7, ncols=8) 
extent(r)<-c(-60,60,-50,50) #avoid touching cells at the west and south edges
r[]<-0
r[4,5]<-1 #value at the upper edge
r[1,3]<-1 #value at the left edge
r[5,1]<-1 #value at the center
plot(r) 

the weight matrix must have a zero value for the own value:

nb.w<-matrix(c(1,1,1,1,0,1,1,1,1),ncol=3)

To take the avereage of all neighboring cells (withoug the own cell value) you can create your own funciton:

mean.B.style<-function(x){sum(x,na.rm=T)/(ncell(nb.w)-1)}
# sum of all values devided by 8 (nr. of neighbors) 
# B style, in reference to the spdep::nb2listw function

To take into account, that edges have fewer neighbors you can adjust the weights with:

mean.W.style<-function(x){sum(x,na.rm=T)/(length(x[!is.na(x)])-1)}
# W style, in reference to the spdep::nb2listw function

With either of these functions you can now create a new raste containing the spatial lags:

nb.r<-focal(r,nb.w,pad=T,NAonly=F,fun=mean.B.style)
plot(nb.r)

base raster plot

spatial lag raster plot