Skip to content

raster: extract() based in a pixel value

2 messages · ASANTOS, Rafael Wüest

#
Dear R-sig-geo Members,

I'd like to use de extract() function using the raster package for 
calculate a proportion of pixel with "1"s values inside a buffer given 
some coordinates in a raster. I try to create a function for this 
without success, in my hypothetical example:

#Package
library(raster)

# Create a raster
ras <- raster(ncol=1000, nrow=1000)
set.seed(0)
values(ras) <- runif(ncell(ras))
values(ras)[values(ras) > 0.5] = 1
values(ras)[values(ras) < 0.5] = NA

# Create some coordinates
pts<-sampleRandom(ras, size=30, xy=TRUE)
pts.df<-as.data.frame(pts)
pts.df$area<-rnorm(30, mean=10)## Here just for create a artificial 
covariate without any direct implication in my question

#Function for extract proportion of 1s values
percentual_1s<- function(x,...) {
 ? leng1<-length(values(x) ==1) # length of 1s pixels values
 ? lengtotal<-length(x) # total length of pixels inside buffer
 ? perc<-(leng1/lengtotal)*100
 ? return(perc)
}

# Extract the desirable proportion in a circular 100000 units buffer
cent_max <- extract(ras, # raster layer
 ??? cbind(pts.df$x,pts.df$y),??? # SPDF with centroids for buffer
 ??? buffer = 100000,???????????? # buffer size
 ??? fun=percentual_1s,?????????? # what to value to extract
 ??? df=TRUE)


Here doesn't work, despite the code look like Ok. My perfect output is:

#??????? x????? y layer????? area?? percentual_1s
#1 -109.26 -43.65???? 1 10.349010?? 23.15
#2?? 93.42 -87.21???? 1? 9.861920?? 45.18
#3?? 57.06? 86.85???? 1? 8.642071?? 74.32
#4 -109.98 -45.63???? 1 10.376485?? 11.56
#5? -92.34? 37.89???? 1 10.375138?? 56.89
#6?? 19.62? 21.51???? 1? 8.963949?? 88.15


Please any ideas or any another package that help in this operation?

Thanks in advanced,
#
Hi there

adapt your function for extraction as follows:

percentual_1s<- function(x,...) {
  leng1<-length(which(x==1))
  lengtotal<-length(x) 
  perc<-(leng1/lengtotal)*100
  return(perc)
}

And add "na.rm = FALSE" in the extract call:

cent_max <- extract(ras, pts[,1:2], buffer = 1000000, fun=percentual_1s, df=TRUE, na.rm = FALSE)

That gives me:

head(cent_max)
  ID    layer
1  1 49.81865
2  2 50.27545
3  3 50.03113
4  4 50.29819
5  5 50.39391
6  6 48.89556

The values I get make sense to me (around 50%).

HTH, Rafi