Search code examples
rrasterr-sfterra

Extract raster values using circular buffer and a custom function

I'd like to calculate the percentual 1´s presence in a binomial raster using circular area using a sf object (buffer arond points). I have tried to use the extract function but I am not able to create a final correct object. I try to:

# Packages
library(sf)
library(raster)
library(ggplot2)

# Create some points
set.seed(1)
df <- data.frame(
  gr = c(rep("a",5),rep("b",5)),
  x  = rnorm(10),
  y = rnorm(10)
)
df <- st_as_sf(df,coords = c("x","y"),remove = F, crs = 4326)
df.laea = st_transform(
  df,
  crs = "+proj=laea +x_0=4600000 +y_0=4600000 +lon_0=0.13 +lat_0=0.24 +datum=WGS84 +units=m"
)

# Create abinomial values raster
r <- raster(extent(df.laea) + 3, resolution=1000)
values(r) <- rep(c(1,0),times=length(r)/2)
r.df <- as.data.frame(r, xy = TRUE)

# Buffer around the points
sample.area <- st_buffer(df.laea,dist=10000)
sample.area$area <- rnorm(n=nrow(sample.area),mean=22)

ggplot() + 
  geom_raster(data = r.df, 
               aes(x = x, y = y, fill = layer)) +
  geom_sf(data = sample.area, fill = "white", color = "white") +
  geom_sf(data = df.laea,color = "red")

test

#Function for extract percentual 1s 
perc1<- function(x,...) {
  leng1<-length(which(x==1)) # length of 1s pixels
  lengtotal<-length(na.omit(x)) # total length of pixels inside each circular sample
  perc<-(leng1/lengtotal)
  return(perc)
}

res <- NULL
for (s in 1: length(sample.area)){
# Extract inside ech circular buffer 
square_cc <- raster::extract(r,             # raster layer
    sample.area[s],   # spatial polygon for extraction
    fun=perc1,         # what to value to extract
    df=TRUE, na.rm = FALSE)         # return a dataframe
res <- rbind(res,c(s,square_cc))         
}
#
res
#      ID         layer     
# [1,] integer,10 numeric,10
# [2,] integer,10 numeric,10
# [3,] integer,10 numeric,10
# [4,] integer,10 numeric,10
# [5,] integer,10 numeric,10

And I'd like as output something like this and if have any NA's (outside the area) I need to remove it:

s  area      square_cc   
1  22.91898  0.56230
...
6 -0.3053884 0.59390132 0.73221

In my example, I loose 3 circular areas a cause of NA's or circle outside the raster. Please, any help with it?

Solution

  • Since you only have 0 and 1 values, you can compute the mean to get the fraction of cells that is 1 (when ignoring NAs).

    You can do that like this:

    library(terra)

set.seed(1)
v <- vect(cbind(x= rnorm(10), y = rnorm(10)), crs="+proj=longlat")
v <- project(v, "+proj=laea +x_0=4600000 +y_0=4600000 +lon_0=0.13 +lat_0=0.24 +datum=WGS84 +units=m")

r <- rast(ext(v) + 3, res=1000, crs=crs(v), vals=c(1,0))

b <- buffer(v, 10000)
e <- extract(r, b, mean, na.rm=TRUE)

# if the buffer is small relative to the raster cells size
e <- extract(r, b, mean, na.rm=TRUE, exact=TRUE)

    As expected for this raster, the values are all close to 0.5

    head(e)
#  ID     lyr.1
#1  1 0.5028999
#2  2 0.5000373
#3  3 0.5020881
#4  4 0.5036896
#5  5 0.5012382
#6  6 0.4996306