无法在 R 中有趣地计算光栅堆栈上的逐像素回归

Question

我正在使用栅格，我有一个 7n 层的 RasterStack。我想使用下面的公式计算逐像素回归。我试图用raster::calc 来做这件事，但是我的函数失败了，并显示了消息：

'lm.fit中的错误（x，y，偏移量=偏移量，singular.ok=singular.ok， ...) : 0 (non-NA) 案例。'

但是所有栅格都可以，并且包含数字（不仅是 NA），我可以绘制它， 我可以用公式计算一般线性回归

 cr.sig=lm (raster::as.array(MK_trend.EVI.sig_Only) ~ raster::as.array(stack.pet)+raster::as.array(stack.tmp)+raster::as.array(stack.vap)+raster::as.array(stack.pre)+raster::as.array(stack.wet)+raster::as.array(stack.dtr))

但是当我用

堆叠图层时

allData = stack(MK_trend.EVI.sig_Only,stack.dtr,stack.wet,stack.pre,stack.vap,stack.tmp,stack.pet)

并尝试计算函数

    # Regression Function, R2
lmFun=function(x){
    x1=as.vector(x);
    if (is.na(x1[1])){
        NA 
    } else {
        m = lm(x1[1] ~ x1[2]+x1[3]+x1[4]+x1[5]+x1[6]+x1[7])
        return(summary(m)$r.squared)
    }
}

我看到了错误消息。
我在 R 和编程方面很新，所以，也许，有一些愚蠢的错误？ 为了使处理工作，我将不胜感激。

Answer 1

您可以使用calc 进行逐像素（局部）回归，但您的公式似乎表明您需要其他东西（全局模型）。

如果回归是像素级的，那么每个单元格的 x 和 y 值的数量将相等，并且可以使用 calc。有关示例，请参阅?calc。

相反，每个单元格有 1 个 y（独立）和 6 个 x（因）变量值。这表明您需要一个全局模型。为此，您可以执行以下操作：

library(raster)
# example data
r <- raster(nrow=10, ncol=10)
set.seed(0) 
s <- stack(lapply(1:7, function(i) setValues(r, rnorm(ncell(r), i, 3))))  
x <- values(s)

# model
m <- lm(layer.1 ~ ., data=data.frame(x))

# prediction
p <- predict(s, m)

这需要将所有值加载到内存中。如果你不能这样做，你可以采取一个大的常规样本。见sampleRegular

并说明为什么您的方法不起作用：

testFun=function(x1){
    lm(x1[1] ~ x1[2]+x1[3]+x1[4]+x1[5]+x1[6]+x1[7])
}

# first cell
v <- s[1]
v
#      layer.1  layer.2   layer.3  layer.4  layer.5  layer.6  layer.7
#[1,] 4.788863 4.345578 -0.137153 3.626695 3.829971 4.120895 1.936597

m <- testFun(v)
m
#Call:
#lm(formula = x1[1] ~ x1[2] + x1[3] + x1[4] + x1[5] + x1[6] + x1[7])

#Coefficients:
#(Intercept)        x1[2]        x1[3]        x1[4]        x1[5]        x1[6]        x1[7]  
#      4.789           NA           NA           NA           NA           NA           NA  


summary(m)$r.squared
# 0

即使我没有收到您报告的错误消息（但所有 R^2 值都为零）。