从逐像素回归中提取残差答案

【问题标题】：extracting residuals from pixel by pixel regression从逐像素回归中提取残差
【发布时间】：2017-10-23 21:04:00
【问题描述】：

我正在尝试在 NDVI/降水的栅格堆栈上逐个像素地从回归运行中提取残差。当我用我的一小部分数据运行它时，我的脚本就可以工作。但是当我尝试运行整个研究区域时，我得到：“setValues(out, x) 中的错误：值必须是数字、整数、逻辑或因子”

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

【问题讨论】：

哼，也许你的公式参数有问题？我会明确地调用列。
指向您数据的链接已过时。

标签： r raster lm

【解决方案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/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)

【讨论】：