r语言 - 输出随机森林回归残差作为单个光栅?

Question

我的目标是使用Kriging方法降低回归模型的残差。

我在R中使用随机森林包来执行回归。我的数据集由一个因变量和一个自变量组成(以栅格格式)。当我导出随机森林(RF)回归的残差时，输出是一个堆叠的栅格(500栅格，作为树的数量)。

我的问题是，我可以为残差导出单个光栅，就像我执行简单线性回归时一样?

我正在使用的代码和数据集:

library(randomForest)
library(terra)
library(raster)
wd = "path/"
ntl = rast(paste0(wd, "ntl.tif"))
covar = rast(paste0(wd, "agbh.tif"))
covar = resample(covar, ntl, method = "bilinear")
s = c(ntl, covar)
names(s) = c("ntl", "covar")
m <- randomForest(ntl ~ ., 
data = as.data.frame(s), 
mtry = 1, 
importance = TRUE, 
na.action = na.omit)
# extract residuals
rsds = s$ntl - predict(m)

数据集:

ntl = rast(ncols=113, nrows=56, nlyrs=1, xmin=4498800, xmax=4544000, ymin=4292800, ymax=4315200, names=c('ntl'), crs='PROJCRS["World_Mollweide",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Mollweide"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]')
covar = rast(ncols=452, nrows=225, nlyrs=1, xmin=4498900, xmax=4544100, ymin=4292600, ymax=4315100, names=c('agbh'), crs='PROJCRS["World_Mollweide",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Mollweide"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]')

Answer 1

示例数据

library(terra)
s <- rast(system.file("ex/logo.tif", package="terra")) [[1:2]] 
names(s) = c("ntl", "covar")

解决方案

library(randomForest)
m <- randomForest(ntl~., data=s, mtry=1, importance=TRUE, na.action=na.omit)
p <- predict(s, m)
rsds <- s$ntl - p
rsds
#class       : SpatRaster 
#dimensions  : 77, 101, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 
#source(s)   : memory
#name        :       ntl 
#min value   : -15.20889 
#max value   :  14.35584