1

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.)

raster artifacts

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!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.