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

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

因为我可以提取斜率和拦截。我只是不能提取残差。

任何想法如何解决这个问题？

这里是我的脚本：

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

来源

2017-10-17 Neija

嗯，也许有一些错误的公式参数？我会明确地调用列。 –

您的数据链接已过时。 – Borealis

你的功能仅在第一层显示NA值，以避免拟合模型试验。但其他层可能有NA。你知道，因为你在lm fit中加了na.action = na.exclude。
问题是，如果模型因为NA而删除了一些值，则残差只会具有非NA值的长度。这意味着您的结果r矢量将具有不同的长度，具体取决于图层中NA值的数量。然后，calc不能够在堆栈中将不同长度的结果组合成定义数量的层。
为了避免这种情况，您需要在函数中指定r的长度，并仅将属性残差指定为非NA值。
我建议以下函数可以在您提供的数据集上工作。我增加了（1）如果你想扩展你的探索（如果你想扩展你的探索），可以比较更多的层次（用nlayers），（2）如果每层只有两个值进行比较（完美模型）如果由于某种原因该模型可以适用，则添加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)

来源

2017-10-23 13:04:00

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

回答

相关问题