R filter raster using focal() with threshold - defining correct function

Question

I have the issue of defining the right function for my task. For filtering an image I have set up a filter matrix eucdis with

library(rgdal)
library(raster)
refm=matrix(1,nrow=11,ncol=11)
M = dim(refm)[1]

N = dim(refm)[2]

eucdis = matrix(NaN, nrow=11, ncol=11)
for (i in -5:5){
      for (j in -5:5){
            eucdis[i+6,j+6] = 2*(sqrt(sum(abs(0-i)^2+abs(0-j)^2))) #euclidean distance of the moving matrix
            eucdis[6,6]=1
            eucdis[eucdis>10]=0
            eucdis[eucdis>0]=1
      }
}

Using the example raster

f <- system.file("external/test.grd", package="raster")
f
r <- raster(f)

I want to filter all values of that raster that have a certain value, say 200 within 10% (=8) of the moving eucdis filter matrix

s=focal(x=r,w=eucdis,fun=function(w) {if (length(w[w==1])>=8) {s=1} else {s=0}})

But this only gives me all values where the eucdis filter matrix has at least 8 pixel with any values of r. If I add the constraint about r[r>=200] it is not working as I thought it would. It is not taking the second constraint into account.

s=focal(x=r,w=eucdis,fun=function(w,x) {
           if (length(w[w==1])>=8 | x[x>=200]){s=1} else {s=0}}) 
# I also tried & and &&

If anyone can help me please. I have spend days already and can't figure it out my self.

Thank you,

Anne

Answer 1

The function passed to focal doesn't refer to the weights matrix. Rather, it refers to the cells of r that fall within the moving window (and these cells' relative contribution to the function's return value is controlled by the weights matrix). So, where you have used function(w) {if (length(w[w==1])>=8) 1 else 0}, you're actually saying that you want to return 1 if the focal subset of r has at least 8 cells with value equal to 1 (and return 0 otherwise).

One to achieve your aim is to perform a focal sum on a binary raster that has been thresholded at 200. The function that you would apply to the moving window would be sum, and the output of this focal sum would indicate the number of cells of your thresholded raster that have value 1 (this corresponds to the number of cells of r that have value >= 200 within the moving window).

library(raster)
r <- raster(system.file("external/test.grd", package="raster"))

m <- matrix(2 * pointDistance(expand.grid(-5:5, -5:5), c(0, 0), lonlat=FALSE),
            ncol=11, nrow=11)
m <- m <= 10

r2 <- focal(r >= 200, m, sum, na.rm=TRUE, pad=TRUE)
plot(r2)

You can then check which cells of that raster have value >= 8.

r3 <- r2 >= 8
plot(r3)

In this case, almost all the cells meet your criteria.