extracting residuals from pixel by pixel regression

Question

I am trying to extract the residuals from a regression run pixel by pixel on a raster stack of NDVI/precipitation. My script works when i run it with a small part of my data. But when i try to run the whole of my study area i get: "Error in setValues(out, x) : values must be numeric, integer, logical or factor"

The lm works, since I can extract both slope and intercept. I just cant extract the residuals.

Any idea of how this could be fixed?

Here is my script:

setwd("F:/working folder/test")
gimms <- list.files(pattern="*ndvi.tif")
ndvi <- stack(gimms)
precip <- list.files(pattern="*pre.tif")
pre <- stack(precip)
s <- stack(ndvi,pre)

residualfun = function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[1:6] ~ x[7:12], na.action=na.exclude)
r <- residuals.lm(m)
return (r)}}

res <- calc(s,residualfun)

And here is my data: https://1drv.ms/u/s!AhwCgWqhyyDclJRjhh6GtentxFOKwQ

Answer 1

Your function only test if the first layer shows NA values to avoid fitting the model. But there may be NA in other layers. You know that because you added na.action = na.exclude in your lm fit.
The problem is that if the model removes some values because of NA, the residuals will only have the length of the non-NA values. This means that your resulting r vector will have different lengths depending on the amount of NA values in layers. Then, calc is not be able to combine results of different lengths in a stack a a defined number of layers.
To avoid that, you need to specify the length of r in your function and attribute residuals only to non-NA values.
I propose the following function that now works on the dataset your provided. I added (1) the possibility compare more layers of each if you want to extend your exploration (with nlayers), (2) avoid fitting the model if there are only two values to compare in each layer (perfect model), (3) added a try if for any reason the model can fit, this will output values of -1e32 easily findable for further testing.

library(raster)
setwd("/mnt/Data/Stackoverflow/test")
gimms <- list.files(pattern="*ndvi.tif")
ndvi <- stack(gimms)
precip <- list.files(pattern="*pre.tif")
pre <- stack(precip)
s <- stack(ndvi,pre)

# Number of layers of each
nlayers <- 6

residualfun <- function(x) {
  r <- rep(NA, nlayers)
  obs <- x[1:nlayers]
  cov <- x[nlayers + 1:nlayers]
    # Remove NA values before model
    x.nona <- which(!is.na(obs) & !is.na(cov))
    # If more than 2 points proceed to lm
    if (length(x.nona) > 2) {
      m <- NA
      try(m <- lm(obs[x.nona] ~ cov[x.nona]))
      # If model worked, calculate residuals
      if (is(m)[1] == "lm") {
        r[x.nona] <- residuals.lm(m)
      } else {
        # alternate value to find where model did not work
        r[x.nona] <- -1e32
      }
    }
return(r)
}

res <- calc(s, residualfun)