从逐像素回归中提取残差
extracting residuals from pixel by pixel regression
我正在尝试在 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)
您的函数仅测试第一层是否显示 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)
我正在尝试在 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)
您的函数仅测试第一层是否显示 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)