我想知道是否可以将栅格包的叠加函数与调用向量列表的函数一起使用，以基于两个栅格执行一些计算。到目前为止，我只看到了函数在不调用外部数据的情况下执行一些栅格代数的示例。
下面，我提供了一些玩具代码来说明我正在尝试做什么，但也可以提供一些关于我真正问题的上下文。具体来说，我需要将每个像素分类为零(不存在(或一个(存在(外壳。住房存在的可能性与覆盖像素的建成区百分比(下面的栅格"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]))})
错误是
不能使用此公式，可能是因为它没有矢量化
除了上面的玩具代码之外，我想单独处理每个分类类并使用函数 calc 而不是函数覆盖(见下文(。但是，对于大型栅格来说，这是非常慢的(如果不是不可能的话(，所以我认为叠加会更好。
# 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)
我最近遇到了同样的错误。我的解决方案是"矢量化"传递给覆盖的函数。这会使用 mapply 创建函数的包装器，以便覆盖可以使用它。我能够通过使用矢量化(见下文(让您的代码运行。
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)
## edits below
# define function
f <- function(x,y,...){return(rbinom(n=1, size=1, prob=probs[[y]][x]))}
# run function using overlay with Vectorize
o <- overlay(r1,r2,fun=Vectorize(f))
这将生成概率的栅格图层。它还产生了以下错误：
在 rbinom(n = 1， size = 1， prob = probs[[y]][x]( 中：生成的 NA
我不确定此错误是否会对您的真实数据造成问题。
您还可以参考此处的第二个答案以获取另一个工作示例。