从逐像素回归中提取残差

Question

我正在尝试在 NDVI/precipitation 的光栅堆栈上逐个像素地从回归运行中提取残差。当我运行它包含我的一小部分数据时，我的脚本可以工作。但是当我尝试运行我的整个学习区域时，我得到： "Error in setValues(out, x) : values must be numeric, integer, logical or factor"

lm 有效，因为我可以提取斜率和截距。我只是无法提取残差。

知道如何解决这个问题吗？

这是我的脚本：

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)

这是我的数据：https://1drv.ms/u/s!AhwCgWqhyyDclJRjhh6GtentxFOKwQ

Answer 1

您的函数仅测试第一层是否显示 NA 值以避免拟合模型。但是其他层可能有NA。您知道这一点，因为您在 lm 中添加了 na.action = na.exclude。
问题是，如果模型因为 NA 而删除了一些值，残差将只有非 NA 值的长度。这意味着您生成的 r 向量将具有不同的长度，具体取决于层中 NA 值的数量。然后，calc 无法将堆栈中不同长度的结果组合到定义的层数。
为避免这种情况，您需要在函数中指定 r 的长度，并将残差仅指定为非 NA 值。
我建议以下功能现在适用于您提供的数据集。我添加了 (1) 如果你想扩展你的探索（使用 nlayers），可以比较每一层的更多层，(2) 如果每一层只有两个值要比较，则避免拟合模型（完美模型） , (3) 添加了一个 try 如果出于任何原因模型可以拟合，这将输出 -1e32 的值以便于进一步测试。

library(raster)
setwd("/mnt/Data/Whosebug/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)

从逐像素回归中提取残差

extracting residuals from pixel by pixel regression

r

raster

lm