带有调用矢量列表的函数的raster::overlay

我想知道是否可以将光栅包的overlay函数与调用矢量列表以基于两个光栅执行某些计算的函数一起使用。到目前为止,我只看到一些函数在不调用外部数据的情况下执行一些光栅代数的示例。

下面我将提供一些玩具代码来说明我想要做什么,但也可以提供一些关于我真正问题的上下文。具体来说,我需要将每个像素分为零(不存在)或一(存在)个外壳。房屋存在的可能性与覆盖像素的建筑面积百分比(下面的光栅r1)和土地覆盖类型(下面的光栅r2)有关。这种可能性是基于参考数据而知道的,参考数据存储在下面类似probs的列表中。

library(raster)

# continuous and categorical maps
r1<-r2<-raster()
r1[]<-round(runif(ncell(r1))*100)
r2[]<-1
r2[1:30000]<-2

# probability of housing presence in each stratum
prob1<-1:100/100 
prob2<-log(1:100)/max(log(1:100))

# list of probabilities to be used in overlay
probs<-list(prob1,prob2)

# overlay - not working
o<-overlay(r1,r2,fun=function(x,y,...){return(rbinom(n=1, size=1, prob=probs[[y]][x]))})

错误是

无法使用此公式,可能是因为它未矢量化

除了上面的toy代码之外,我还想分别处理每个分类类,并使用函数calc而不是函数overlay(见下文)。不过,这是非常缓慢(如果不是不可能)的大光栅,所以我认为覆盖会更好。

# alternative: loop across categorical classes (extremely slow for large rasters)
r<-list()
for(i in 1:2){
  stratum<-r2
  stratum[Which(stratum !=i)]<-NA
  r[[i]]<-calc(r1, fun=function(x,...){return(rbinom(n=1, size=1, prob=probs[[i]][x]))})
  r[[i]]<-mask(r[[i]],stratum)
}

r<-stack(r)
r<-sum(r,na.rm=T)

par(mfrow=c(1,3))
plot(r1)
plot(r2)
plot(r)