Rolling (moving) window linear regression - extract residuals

Question

I want to perform linear regression with a rolling (moving) window with size 5 using raster data.

I tried a code using the raster package and the function localFun, and when I am trying to export the residuals I am getting the following error: Error in setValues(x, value) : values must be numeric, logical or factor. Here is the code:

library(raster)

ntl = raster("path/ntl.tif") # dependent variable
tirs = raster("path/tirs.tif") # independent variable
tirs_aggr = resample(tirs, ntl, method = 'bilinear')

rfun1 <- function(x, y, ...) {
  d <- na.omit(data.frame(x, y))
  if (nrow(d) < 5) return(NA)
  m <- lm(y~x, data = d)
  # return intercept
  coefficients(m)[1]
}

rfun2 <- function(x, y, ...) {
  d <- na.omit(data.frame(x, y))
  if (nrow(d) < 5) return(NA)
  m <- lm(y~x, data = d)
  # return slope
  coefficients(m)[2]
}

rfun3 <- function(x, y, ...) {
  d <- na.omit(data.frame(x, y))
  if (nrow(d) < 5) return(NA)
  m <- lm(y~x, data = d)$residuals # doesn't work
  # return residuals
  # residuals(m) or m$residuals # doesn't work
}

ff = localFun(ntl, tirs_aggr, ngb = 5, fun = rfun1)
ff2 = localFun(ntl, tirs_aggr, ngb = 5, fun = rfun2)
ff3 = localFun(ntl, tirs_aggr, ngb = 5, fun = rfun3) # returns error

My question is this: How can I extract the coefficients and the residuals as raster layers from a moving window linear regression?

The data:

ntl = raster(ncols=116, nrows=98, xmn=509587.9392, xmx=550187.9392, ymn=161637.6238, ymx=195937.6238, crs='+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs')

tirs = raster(ncols=409, nrows=344, xmn=509600, xmx=550500, ymn=161700, ymx=196100, crs='+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs')

Answer 1

Check package's code. It's not possible to return a raster with 2 or more layers:

out <- raster(x)

Since your matrix is 5x5, you can get up to 25 values from m$residuals or residuales(m), which is incompatible with raster package.

You can do it not directly with terra. I give you an example:

library(terra)

# example rasters
ntl <- rast(nrows=18, ncols=36)
tirs_aggr <- ntl

set.seed(123)
values(ntl) <- rnorm(n=ncell(ntl))
values(tirs_aggr) <- rnorm(n=ncell(ntl))

#get focal values for each window
fntl <- focalValues(ntl, w=5, row=1, nrows=nrow(ntl), fill=NA)
ftirs_aggr <- focalValues(tirs_aggr, w=5, row=1, nrows=nrow(ntl), fill=NA)

# funtion to compute all you want
rfun_total <- function(x) {
  x1 <- x[1:(length(x)/2)]
  y1 <- x[((length(x)/2)+1):length(x)]
  d <- na.omit(data.frame(x1, y1))
  
  result <- rep(NA, length.out = 27)
  names(result) <- c('inter','x',paste0('res',1:25))
  
  if (nrow(d) >= 5){
    m <- lm(y1~x1, data = d)
    # return slope
    result[1:(length(m$residuals)+2)] <- c(coefficients(m),m$residuals)
    
  }
  
  result
}

# concatenate results and apply function
result <- apply(cbind(fntl,ftirs_aggr), MARGIN = 1, FUN = rfun_total)

final_raster <- rast(ntl, nlyrs=27)

rlist <- apply(result, 1, function(x) setValues(ntl,x))

rast(rlist)

# class       : SpatRaster 
# dimensions  : 18, 36, 27  (nrow, ncol, nlyr)
# resolution  : 10, 10  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 
# source(s)   : memory
# names       :      inter,          x,      res1,      res2,      res3,      res4, ... 
# min values  : -0.5300309, -0.8214684, -2.699939, -2.736957, -2.534735, -2.648687, ... 
# max values  :  0.7583305,  0.6152796,  3.283922,  3.099785,  3.067819,  3.044425, ... 

Also, you can add your need to this thread. I find you placed a good point here

Finally, with this example you computed the linear model once, no need to do it so many times as you have been doing so far