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,   
type='circle')

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
}
return(Ratio)
}

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

library(raster)
library(sp)

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
fw

( rmean <- focal(r, w = fw ,fun = mean, pad=TRUE, 
                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 <- focalWeight(r, 810, type='circle')
fw[fw > 0] <- 1
( rratio <- focal(r, w = fw ,fun = ratio, pad=TRUE, 
                padValue = 0) )
  • 1
    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 at 18:21
  • @Spacedman, thanks, i forgot that they are weights that sum to 1 and not a Gaussian distribution. – Jeffrey Evans Aug 10 at 20: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.

library(raster)
#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

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.