0

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)

#read in raster
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.

2
  • A percentage would be in range 0-100.
    – Vince
    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
    – seak23
    Commented Aug 13, 2023 at 3:24

1 Answer 1

1

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
1
  • 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?
    – seak23
    Commented Aug 14, 2023 at 5:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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