I want to calculate the ratio of two different classes in a raster with a circular moving window. My classes are 0 and 1 (I do not want to include the NAs). The radius should be 5 meter. I thought my code would work in R until i tried it out in Esri with focal statistics in order to compare the results and got completely different (but a lot more plausible) results. What might be the problem with my code in R? As far as I can see they should be doing the same thing. Unfortunatley I could not figure out the exact code that Esri uses...

My raster contains only integers of 1, 0, NA

fw <- focalWeight(modell, 5,   

result <- focal(modell, w = fw ,fun = 
mean, pad=TRUE, padValue = 0, na.rm = TRUE)

I was not 100% sure how the mean-function would handle 0s and NAs, so I also tried to write a function myself and change all the NAs in my code to -1 . It did not improve the output though.

modell[is.na(modell[])] <- -1

fw <- focalWeight(modell, 5, type='circle')

RatioFun <- function(fw){
  NCells <- sum(fw >= 0)
  SumCells <- sum(fw[which(fw >= 0)])
  if (NCells == 0){
    SumCells <- 0
  } else {
    Ratio <- SumCells/NCells

focal(modell, w = fw ,fun = RatioFun,  
pad=TRUE, padValue = 0, na.rm = TRUE)

This is how i use the focal statistics in ArcGis. enter image description here

This are the two different results. enter image description here

  • Interesting. I reproduced differences with success... Something is not working in R – aldo_tapia Aug 10 '18 at 13:32

The result of the raster::focalWeight is a Gaussian kernel, thus the reference to "weights" in the function name. If you simply coerce it to a binary matrix then you should get the expected results.

Let's create some example data


r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
            ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
r[] <- rpois(ncell(r), lambda=1)
r <- calc(r, fun=function(x) { x[x >= 1] <- 1; return(x) } )

Here we create the focalWeight object for 810 meters, and take a look at the resulting object (remember that the d argument is in projection distance units and not number of cells as in a standard rectangular window). If we apply it to the focal function then we get unexpected results.

( fw <- focalWeight(r, 810, type='circle') )

( rmean <- focal(r, w = fw ,fun = mean, pad=TRUE, 
                padValue = 0, na.rm = TRUE) )

Now, let's coerce the focal matrix so it is binary and apply it to the raster. Presto, the results now make sense. The raster values associated with the zero values in the focal matrix are not passed to the function called within raster::focal.

fw[fw > 0] <- 1

                padValue = 0, na.rm = TRUE) )

As a gut check, we can also write our own ratio function and pass it to focal.

ratio <- function(x, value = 1) {
  x <- x[!is.na(x)]
    n1 <- length(x[x == value])
  if( length(x) == 0 ) {
    p <- NA
  } else if( length(length(x[x != value])) == 0 ) {
      p = 1   
    } else if( length(n1) == 0 ) {
      p = 0
    } else {
      p <- n1 / length(x)
  return ( p )

fw[fw > 0] <- 1
( rratio <- focal(r, w = fw ,fun = ratio, pad=TRUE, 
                padValue = 0) )
  • 2
    focalWeight with type="circle" isn't a Gaussian kernel, but it does integrate to 1, so the values are all constant and 1/number of non-zero cells – Spacedman Aug 10 '18 at 18:21
  • 1
    @Spacedman, thanks, i forgot that they are weights that sum to 1 and not a Gaussian distribution. – Jeffrey Evans Aug 10 '18 at 20:15
  • So changing the matrix does the trick. I should have thought about it myself, but somehow did not make that connection... With focalWeight and type='circle' the cells can be 0 or sum up to 0 (which means they are often really small numbers). After changing them to only cells of 0 and 1 it works perfectly. fw[fw > 0] <- 1 Thanks again and sorry for the last confirmation! – Li12 Aug 23 '18 at 14:15

I have had similar problems myself. I used the sum function instead and it lines up with ArcMap's focal mean, +/- some decimal rounding.

#Binary raster, 1 = wildfire, 0 = no wildfire, cellsize = 30m, NAD83 UTM11N
wf <- raster("H:/testingFocal/narWildfire.tif")
#Generate 500m moving window
fw <- focalWeight(wf,500,type = "circle")
#Execute focal sum
wfFocalSum <- focal(wf, fw, fun=sum, na.rm=FALSE)
#Write focal sum to disk to compare against ArcMap focal mean
writeRaster(wfFocalSum, "H:/testingFocal/R500m.tif", overwrite = TRUE)

Execute ArcMap focal mean statistics on same binary wildfire raster:

enter image description here

Check both outputs with the identity tool in ArcMap to verify that values are the same.

enter image description here

  • Thanks for the tipp. I see your rounding might be different. The values also do not differ so much. My values were really different. I was also thinking if some rounding might change the outcome, but it could not have influenced my results that much. Thanks for the hint though. It is a good point to check! – Li12 Aug 23 '18 at 13:55

Your Answer

