How can I get the correct mean from focal() with a circular window?

Question

I am trying to perform focal statistics on a raster in R using a circular window (defined with focalWeight) and the focal function. In my specific case, I want each cell to be the mean of all cell values within a 100 m radius. Because I have NAs in my dataset, I followed the advice of Robert Hijmans in this post, which is to give weights a value of 1 and use fun = mean rather than the default sum function.

This seemed a perfect solution; however, the values produced are wrong. (This post describes a similar problem but unfortunately never received an answer.)

After much reading and testing I have come to realise that the problem is linked to the denominator used for the calculation of the mean. My question is, can I achieve the correct result without resorting to a clunky, slower custom function?

Here is a small example of the problem (demonstrated without NAs to keep it simple). Let's create a raster where all cells have a value of 1.

# Create a raster of 9 x 9 cells (res 1 m when projected in GB national grid)
r<-raster(extent(0,9,0,9), ncol = 9, nrow = 9, vals = 1, crs = "+init=epsg:27700")   # all the cells have value of 1

# Set circular window of 2m radius
fw <- focalWeight(r, 2, "circle")

Because there are no NA values, we can use the sum function and the weights. As expected, all cells (except edges) receive the correct mean value of 1.

focal_r <- focal(r, fw)

summary(values(focal_r))  # all values of 1

I would expect that changing the weights to values of 0 (outside circle) and 1 (inside circle) and using the mean function would yield the same result (as suggested in R. Hijmans' answer). However:

fw[fw > 0] <- 1   # replacing weights by 1

focal_mean_r <- focal(r, fw, fun = mean) 

summary(values(focal_mean_r))  # now highest value is 0.52

The highest value we can achieve is 0.52, which I found out is the ratio of "1" cells to all cells in the focal window:

ncell(fw[fw > 0])/ncell(fw)

In other words, the function sums the raster values for cells falling in the "1" circular window, but divides the result by the total dimensions of the weights matrix. But surely this is not a desirable behaviour? The 0 cells (and NA cells for that matter, when na.rm is selected) should not count in the denominator. So is this an issue with the method? Or am I doing something wrong?

I even tried to coerce the focal window 0s to NA as suggested in method #4 of this blog post, but it doesn't work and I get a blank raster.

So my workaround is a custom function where I multiply the weights of 1 by the cell values, sum them all, and divide by the number of non-zero cells of the focal window:

custom_r <- focal(r, fw, fun = function(x){sum(x * fw, na.rm = TRUE)/ncell(fw[fw > 0])})

summary(values(custom_r))  # correct results with all values of 1

For the visual folks, you can see artifacts arising from the wrong calculation (left) when I view my raster in QGIS, while the output raster looks normal and smooth using my custom function (right). (NB: My large raster currently has only values of 1 or NA, although values could in practice range from 0 to 2.)

If anyone can confirm that this is the only way to go, I'll live with it. But the custom function runs about 3x slower than the built-in mean method, so if I'm missing something I would really like to know!

Answer 1

I had the same problem using the raster::focal function with a circular mean. A good work-around is to just use the equivalent function in Terra, which is the successor to the Raster package. Here's a continuation of your example (converting to Terra only for the relevant step and then back to the Raster package format, although I'd recommend switching over to Terra completely):

# Create a raster of 9 x 9 cells (res 1 m when projected in GB national grid)
r <- raster(extent(0,9,0,9), ncol = 9, nrow = 9, vals = 1, crs = "+init=epsg:27700")   # all the cells have value of 1

#convert the original Raster package format to the Terra package format 
r_spat_rast <- terra::rast(r)     

# Set circular window of 2m radius 
fw <- focalWeight(r, 2, "circle") #The equivalent Terra function is focalMat
fw[fw > 0] <- 1   # replacing weights by 1

#using the Terra version of the focal function
focal_r_spat_rast <- terra::focal(x = r_spat_rast, w = fw, fun = "mean") 

focal_r <- raster::raster(focal_r_spat_rast ) #converting back to the Raster package format, in case the rest of your workflow is still with the Raster package

summary(values(focal_r))  # all values of 1

Answer 2

I posted this as an issue on github. It is now fixed in the development version of raster (3.5.8) if you use NA's instead of 0's in the focal window. 0 means that the values are multiplied by 0 and that is used in the calculation where as NA indicates to ignore those values. terra::focal does however behave differently than raster::focal when mean is used as the function, as it ignores the 0 weights (i.e. using NA or 0 weights is the same when using mean in terra::focal.

https://github.com/rspatial/raster/issues/218

Answer 3

A somewhat simpler function (but perhaps no faster, though I haven't checked) using raster is to first convert the focal weights matrix to 1s and NAs, then (using focal) multiply the window (x) by these fw values and finally take the mean (removing any NAs in fw or the raster using na.rm=T). This does create NAs at edges of the image with a width of the radius used; padding is necessary (not shown) to deal with those.

fw <- focalWeight(r, 2, "circle")
fw[fw>0] <- 1
fw[fw==0] <- NA 

custom_r <- focal(r, fw, fun = function(x) mean(x * fw, na.rm=T)) 

Warnings are returned but the process seems to work.

EDIT: I just ran the analogous approach using both raster terra and found the same result in a few seconds vs. a couple minutes using the approach I noted above. Critical here is to set na.rm = T.

fw <- focalWeight(r, 2, "circle")
fw[fw>0] <- 1
fw[fw==0] <- NA 

These return the same results (looks like the newest update has fixed the denominator issue, thanks!)

focal_r <- terra::focal(r, fw, fun = "mean", na.rm=T)
focal_r <- raster::focal(r, fw, fun = "mean", na.rm=T)