quantForestError() 对栅格堆栈的预测？

Question

有聪明的人知道如何在栅格堆栈而不是矩阵或 data.frame 上编写 quantForestError() 预测吗？例如...

library(ranger)
library(tidyverse)
library(raster)
library(forestError)
data(iris)
glimpse(iris)

mod<-ranger(Sepal.Length~. -Species,iris,keep.inbag =TRUE,seed=4)#random forest model

#trying to get prediction, and upper and lower prediction interval
#using quantForestError(model,x-training data, new data, y-training data, what I want)$which thing as output
quantForestError(mod,iris[,2:4],iris[10,2:4],iris[,1],what=c("interval"))$pred
#4.823236
quantForestError(mod,iris[,2:4],iris[10,2:4],iris[,1],what=c("interval"))$upper_0.05
#5.091271
quantForestError(mod,iris[,2:4],iris[10,2:4],iris[,1],what=c("interval"))$lower_0.05
#4.293341

#generate raster data
sw <- raster(ncol=10, nrow=10)
set.seed(4)
values(sw) <- runif(ncell(sw),min=iris$Sepal.Width,max=iris$Sepal.Width)
pl <- raster(ncol=10, nrow=10)
set.seed(4)
values(pl) <- runif(ncell(pl),min=iris$Petal.Length,max=iris$Petal.Length)
pw <- raster(ncol=10, nrow=10)
set.seed(4)
values(pw) <- runif(ncell(pw),min=iris$Petal.Width,max=iris$Petal.Width)
all<-stack(sw,pl,pw)
names(all)<-c("Sepal.Width","Petal.Length","Petal.Width")
plot(all)

##try on raster --> error
quantForestError(mod,iris[,2:4],all,iris[,1],what=c("interval"))$pred
#Error in checkXtrainXtest(X.train, X.test) : 
#'X.test' must be a matrix or data.frame of dimension 2


#for what it's worth this is how I would get prediction on raster normally with ranger
plot(predict(all, mod, type='se', seed=4, fun = function(model, ...) predict(model, ...)$predictions))
plot(predict(all, mod, type='se', seed=4, fun = function(model, ...) predict(model, ...)$se))

我确实知道我可以使用护林员预测中的这些预测来获得置信区间，但请参阅问题here。我认为我需要做的是将栅格转换为矩阵，进行预测，然后转换回栅格，但我不确定如何执行此操作。感谢您提供任何意见。

Answer 1

这是将栅格堆栈转换为data.frame，计算预测和预测间隔，然后转换回栅格的方法。它非常慢，但似乎比使用

ranger

包计算标准误差要好。

all2<-as.data.frame(all,xy=TRUE)%>%drop_na()
colnames(all2)
out<-quantForestError(mod,iris[,2:4],all2[,3:5],iris[,1],what=c("interval"))
pred<-cbind(all2[,1:2],out[,1])
lower<-cbind(all2[,1:2],out[,2])
upper<-cbind(all2[,1:2],out[,3])
pred_rst <- rasterFromXYZ(pred,crs=crs(all), res=res(all))
pred_lwr <- rasterFromXYZ(lower,crs=crs(all), res=res(all))
pred_upr <- rasterFromXYZ(upper,crs=crs(all), res=res(all))
plot(stack(pred_lwr,pred_rst,pred_upr))