# Circular focal means with raster and NAs

I am trying to understand circular focal means using rasters having NA values. I made a couple of scenarios to better understand how to combine weight matrices and functions to summarize values. I am using a circular weight matrix considering the focal cell and four neighboring cells.

``````library(raster)

r<-raster(extent(0,1,0,1),ncol=5,nrow=5,vals=1:25)
wm<-focalWeight(r,d=res(r)[1]*1.001,type="circle") # original weight matrix
w<-wm
w[w>0]<-1

w
[,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    1    1
[3,]    0    1    0
``````

Using a first example without NA values and if we focus on the middle value of the raster, (A) I'm not able to correctly compute it using mean and weights of 0/1. It seems that 0 values are included in the calculation as the denominator in the mean calculation is the number of values in the weight matrix (9). The documentation warns about the fact that using weights different from 1 in combination with `na.rm = TRUE` may lead to problems in calculations, but in this case I would have expected to get the correct value. (B) I can compute the correct value by using instead weights of 0/0.2 and the sum function.

``````par(mar=c(1,1,5,1),mfrow=c(2,2))
plot(r,axes=FALSE,legend=FALSE,main="Initial raster")
text(coordinates(r),labels=getValues(r))

### A
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=mean)
plot(r2,axes=FALSE,legend=FALSE,main="A\nMean\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+13+14+18) / 9 = 7.22",cex=0.7)

### B
w<-wm
r2<-focal(r,w=w,fun=sum)
plot(r2,axes=FALSE,legend=FALSE,main="B\nSum\nweights = 0/0.2")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+13+14+18) * 0.2 = 13",cex=0.7)
``````

In this second example, I turn the middle value to NA and try to compute the missing value using the same focal area. As in the first example, the mean method (C) does not work and the denominator is now 8 since the missing value was removed. The sum method (D) does not work anymore, because the weights are not adjusted for the removed NA value. A solution (E) could be to use a 0/1 matrix and my own mean function to remove all 0 weight values before the mean calculation.

``````r[13]<-NA

par(mar=c(1,1,5,1),mfrow=c(2,2))
plot(r,axes=FALSE,legend=FALSE,main="Initial raster with NA")
text(coordinates(r),labels=getValues(r))

### C
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=mean,na.rm=TRUE)
plot(r2,axes=FALSE,legend=FALSE,main="C\nMean with na.rm = TRUE\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) / 8 = 6.5",cex=0.7)

### D
w<-wm
r2<-focal(r,w=w,fun=sum,na.rm=TRUE)
plot(r2,axes=FALSE,legend=FALSE,main="D\nSum with na.rm = TRUE\nweights = 0/0.2")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) * 0.2 = 10.4",cex=0.7)

### E
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=function(i){mean(i[i>0],na.rm=TRUE)})
plot(r2,axes=FALSE,legend=FALSE,main="E\nCustom mean function\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) / 4 = 13",cex=0.7)
``````

The problem with this solution is that (1) it relies on my raster not having values of 0 and (2) using my own function greatly increases computation time. The actual matrix I need to run this on is about 13000 x 13000 pixels with lots of NAs.

``````n<-2000
r<-raster(extent(0,1,0,1),ncol=n,nrow=n,vals=1:(n^2))
w<-wm
w[w>0]<-1
system.time(r2<-focal(r,w=w,fun=mean,na.rm=TRUE))
user  system elapsed
0.29    0.00    0.30
system.time(r2<-focal(r,w=w,fun=function(i){mean(i[i>0],na.rm=TRUE)}))
user  system elapsed
74.42    0.40   76.28
``````

My question is this: I there a more efficient way than using a custom-made function to correctly compute circular focal means when NA values are present?

Checking the documentation with `?focal` function I found that `focalWeight`has the posibility to fill NA values. So you can do just this:

``````library(raster)

r<-raster(extent(0,1,0,1),ncol=5,nrow=5,vals=1:25)
wm<-focalWeight(r,d=res(r)[1]*1.001,type="circle", fillNA = TRUE) # original weight matrix

r[13]<-NA

par(mar=c(1,1,5,1),mfrow=c(2,2))
plot(r,axes=FALSE,legend=FALSE,main="Initial raster with NA")
text(coordinates(r),labels=getValues(r))

### C
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=mean, na.rm=TRUE)
plot(r2,axes=FALSE,legend=FALSE,main="C\nMean with na.rm = TRUE\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) / 8 = 6.5",cex=0.7)
``````

``````# benchmark
n<-2000
r<-raster(extent(0,1,0,1),ncol=n,nrow=n,vals=1:(n^2))
w<-wm
w[w>0]<-1
system.time(r2<-focal(r,w=w,fun=mean,na.rm=TRUE))

> system.time(r2<-focal(r,w=w,fun=mean,na.rm=TRUE))
user  system elapsed
0.065   0.004   0.069
``````

I hope it works for you.