# Calculating statistic for pixels within specified neighborhood in R

I have a raster with coded values of either 0 or 1. For each pixel, I want to calculate the percentage of 1s within a 500m radius. Using the focal() function in the terra package seems to be my best option, but I'm not getting the results I expect.

``````library(terra)

r <- rast("Raster/path")

# Set circular window of 500m radius
fw <- focalMat(r, 500, "circle")
fw[fw > 0] <- 1   # replacing weights by 1

#using focal to get the sum
r_sum <- terra::focal(x = rast(raster_data), w = fw, fun = "sum")

r_prop <- r_sum/7854
``````

This is as far as I got. 7854 is the total number of pixels that can fit in a 500 meter circle. I was expecting the value of each pixel to be a percentage (0-100). Right now I'm only getting a max value of 2, so i think something is off.

• A percentage would be in range 0-100. Commented Aug 12, 2023 at 19:47
• Hey yes sorry I fixed that in my post above. My maximum value is still only around 2 so I'm not quite sure what is going on Commented Aug 13, 2023 at 3:24

The focal weight matrix `fw` is a weight matrix that sums to 1. Here's a tiny example:

``````> m = rast()
> m[]=1
> fw = focalMat(m, 2, "circle")
> fw
[,1]       [,2]       [,3]       [,4]       [,5]
[1,] 0.00000000 0.00000000 0.07692308 0.00000000 0.00000000
[2,] 0.00000000 0.07692308 0.07692308 0.07692308 0.00000000
[3,] 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308
[4,] 0.00000000 0.07692308 0.07692308 0.07692308 0.00000000
[5,] 0.00000000 0.00000000 0.07692308 0.00000000 0.00000000
``````

Then when you do `focal` with that as a weight you get the cell values multiplied by the weight. This is so that the output raster has roughly the same integrated "volume" as the input. For example:

``````> m[]=runif(ncell(m))
> sum(m[])
[1] 32321.5
> ms = focal(m, w=fw, fun="sum")
> sum(ms[],na.rm=TRUE)
[1] 31608.27
``````

If you want the total sum of all cells, then set the values in the weights matrix to 1. As illustration, I'll make a raster of all 2 values, and apply a 3x3 weights matrix - that should return 18 all over:

``````> fw = matrix(1, nrow=3,ncol=3)
> r = rast(matrix(2,6,8))
> focal(r, w=fw, fun="sum")
class       : SpatRaster
dimensions  : 6, 8, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 8, 0, 6  (xmin, xmax, ymin, ymax)
coord. ref. :
source(s)   : memory
name        : focal_sum
min value   :        18
max value   :        18
``````

Great. So for your circle weight matrix, compare with 0 and make it numeric and it should look like this (adjust call to suit your data):

``````> fw = 0+(focalMat(m, 2, "circle") > 0)
> fw
[,1] [,2] [,3] [,4] [,5]
[1,]    0    0    1    0    0
[2,]    0    1    1    1    0
[3,]    1    1    1    1    1
[4,]    0    1    1    1    0
[5,]    0    0    1    0    0
``````
• Thank you. I have attempted this with my raster and got a raster with a min of 0 and maximum of 15690 (similar to the output from ArcGIS Pro's "Focal Statistics"). Now should I divide this area by the area of a circle with 500 meter radius? Commented Aug 14, 2023 at 5:42