library(raster)
library(spdep)
library(rworldmap)
##example country
d<-countriesLow
d<-d[d$ISO3=="IDN",]
# plot(d)
d.indo<-d
##example raster
r<-raster::raster()
extent(r)<-extent(d.indo)
res(r)<-4
#raster values
v<-rep(0,ncell(r))
v[27]<-1
v[29]<-NA
r[]<-v
{plot(r)
plot(d.indo,add=T)}
spdep package#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
v<-lag.listw(nb.w,values(r),zero.policy=F,NAOK=T)
## Warning in lag.listw(nb.w, values(r), zero.policy = F, NAOK = T): NAs in lagged
## values
##new raster
nb.r<-r
nb.r[]<-v
nb.r[]
## [1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [13] 0.000 0.000 0.125 0.125 NA NA NA 0.000 0.000 0.000 0.000 0.000
## [25] 0.000 0.125 0.000 NA 0.000 NA 0.000 0.000 0.000 0.000 0.000 0.000
## [37] 0.200 0.200 NA NA NA 0.000 0.000 0.000
sum(nb.r[],na.rm=T)
## [1] 0.775
{plot(nb.r)
plot(d.indo,add=T)}
The output raster shows higher values at the edges.
nb2listw does not account for the missing values at the edge of a raster.
raster package#neighbor matrix
nb.w<-matrix(c(1,1,1,1,0,1,1,1,1),ncol=3)
#weights function
mean.W.style<-function(x){sum(x,na.rm=T)/(ncell(nb.w)-1)}
#new raster
nb.r<-focal(r,nb.w,pad=T,NAonly=F,fun=mean.W.style)
nb.r[]
## [1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [13] 0.000 0.000 0.125 0.125 0.125 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [25] 0.000 0.125 0.000 0.125 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [37] 0.125 0.125 0.125 0.000 0.000 0.000 0.000 0.000
sum(nb.r[],na.rm=T)
## [1] 1
{plot(nb.r)
plot(d.indo,add=T)}
For more ideas and exmpales see the OSGEO-GRASS application documentation of r.neighbors here.
#neigbor matrix
nb.w<-matrix(ncol=5,nrow=5)
nb.w[]<-0.5
nb.w[3,3]<-0
nb.w[1,1]<-0
nb.w[5,5]<-0
nb.w[5,1]<-0
nb.w[1,5]<-0
nb.w[2,2]<-1
nb.w[2,3]<-1
nb.w[2,4]<-1
nb.w[3,2]<-1
nb.w[3,4]<-1
nb.w[4,2]<-1
nb.w[4,3]<-1
nb.w[4,4]<-1
nb.w[]<-nb.w/sum(nb.w)
sum(nb.w)
## [1] 1
nb.w
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.00000000 0.03571429 0.03571429 0.03571429 0.00000000
## [2,] 0.03571429 0.07142857 0.07142857 0.07142857 0.03571429
## [3,] 0.03571429 0.07142857 0.00000000 0.07142857 0.03571429
## [4,] 0.03571429 0.07142857 0.07142857 0.07142857 0.03571429
## [5,] 0.00000000 0.03571429 0.03571429 0.03571429 0.00000000
#weights function
mean.W.style<-function(x){sum(x,na.rm=T)/(length(nb.w[nb.w>0])-1)}
#new raster
nb.r<-focal(r,nb.w,pad=T,NAonly=F,fun=sum,na.rm=T)
nb.r[]
## [1] 0.00000000 0.00000000 0.00000000 0.03571429 0.03571429 0.03571429
## [7] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [13] 0.00000000 0.03571429 0.07142857 0.07142857 0.07142857 0.03571429
## [19] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [25] 0.03571429 0.07142857 0.00000000 0.07142857 0.03571429 0.00000000
## [31] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.03571429
## [37] 0.07142857 0.07142857 0.07142857 0.03571429 0.00000000 0.00000000
## [43] 0.00000000 0.00000000
sum(nb.r[],na.rm=T)+sum(nb.r[1,])
## [1] 1
{plot(nb.r)
plot(d.indo,add=T)}