Export random forest regression residuals as a single raster?-CodePudding

My goal is to downscale the residuals of a regression model using Kriging methods.

I'm using the randomForest package in R to perform regression. My dataset consists of one dependent and one independent variable (in a raster format). When I export the residuals of the random forest (RF) regression the output is a stacked raster (500 rasters, as the number of trees).

My question is, can I export a single raster for the residuals, like I do when I perform simple linear regression?

The code and the dataset I am using:

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)

The dataset:

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]]]]')

CodePudding user response：

Example data

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

Solution

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