光栅图像的线性回归 - lm 抱怨 NA

Question

我确信这可以用几个字节来解决，但我已经在这个简单的事情上花费了几个小时并且无法摆脱它。我不经常使用 R。

我有 5 个代表 5 个光栅图像的 asciigrid 文件。一些像素确实有值，其他像素确实有 NA。例如，第一张图片可能是这样的：

NA  NA  NA  NA  NA
NA  NA  2   3   NA
NA  0.2 0.3 1   NA
NA  NA  4   NA  NA

第二个可能是：

NA  NA  NA  NA  NA
NA  NA  5   1   NA
NA  0.1 12  12  NA
NA  NA  6   NA  NA

如您所见，NA 位置始终相同，我对此 100% 肯定。我愿意做什么：

• 用read.asciigrid()读取文件；
• 使用raster 包中的values() 在长数组中获取它们的值；
• 创建一个有 5 行的矩阵，每行保存对应映射的值；
• 线性拟合每一列并得到系数。每列代表一个像素，并有 5 个值对应于 5 个地图。
• 使用系数值创建两个新的光栅图像。

我被困在lm。具体来说，它说：Error in lm.fit(...): 0 (non-NA) cases。但是，根据我对输入图的了解，应该有 all NA 的列或根本没有 no NA 的列，如下所示：

NA   NA   NA   NA   0.2  2    NA  ... (lots of other columns)
NA   NA   NA   NA   2    2.1  NA
NA   NA   NA   NA   3    0.5  NA
NA   NA   NA   NA   12   6    NA
NA   NA   NA   NA   0.4  2    NA

我希望输出是：

NA   NA   NA   NA   ..   ..   NA

所以我可以使用系数创建一个新的光栅图像并保持 NA 位置。我哪里错了？在下面粘贴我的代码。谢谢。

library(sp)
library(raster)
library(fields)
names = c('...','...','...','...','...')
x = c(10,20,30,40,50)
x = log(x)
y = vector('list',length=length(x))
rasters = vector('list',length=length(x))
for (name in names) {
  ind = which(name == names)
  rasters[ind] = read.asciigrid(name)
  rasters[ind] = raster(rasters[[ind]])
  y[[ind]] = values(rasters[[ind]])
}

y = t(simplify2array(y))
lModel = lm(y ~ x) // Error here!

这是str(y)的输出：

num [1:5, 1:1260630] NA NA NA NA NA NA NA NA NA NA ...（有时这里会有数字）

编辑

感谢@RobertH，我了解了raster::stack 和raster::calc。我试过了：

x <- log(c(10,20,30,40,50))
fun <- function(y) { lm(y ~ x)$coefficients }
r <- calc(s, fun)

在.calcTest 呼叫中获得一个不起眼的Cannot use this function。我看了raster:::.calcTest 没有成功。我尝试管理所有y 值都是NA 的情况，如下所示：

fun = function(y) { 
  if (any(!is.na(y))) { 
    lm(y ~ x)$coefficients
  } else { 
    NA
  }
}
r <- calc(s,fun)

现在它工作了几分钟，但后来我得到了Error in setValues(out, x) : values must be numeric, integer, logical or factor。但是，通常将 NA 设置为栅格值！我看不出这里有什么问题。

Answer 1

这是获取栅格数据的方法

library(raster)
names = c('...','...','...','...','...')
s <- stack(names)
y <- values(s)

你现在可以做这样的事情了。

x <- log(c(10,20,30,40,50))
# need to exclude the rows that are all NA
i <- rowSums(is.na(y)) < ncol(y)
coef <- apply(y[i, ], 1, function(y) lm(y ~ x)$coefficients)
aa <- matrix(NA, ncol=2, nrow=length(i))
aa[i, ] <- coef
b <- brick(s, nl=2)
values(b) <- aa

但您不需要这样做。要进行这样的回归，我会这样做

fun <- function(y) { lm(y ~ x)$coefficients }
r <- calc(s, fun)

但是因为你的单元格只有 NA 值（跨层），所以这将失败（就像上面的应用一样）。你需要编写一个函数来捕捉这些情况：

funa <- function(y) { 
    if(all(is.na(y))) {
        c(NA, NA)
    } else {
        lm(y ~ x)$coefficients 
    }
}
r <- calc(s, funa)

或者更快的方法

X <- cbind(1, y)
invXtX <- solve(t(X) %*% X) %*% t(X)
quickfun <- function(i) (invXtX %*% i)
m <- calc(s, quickfun) 
names(m) <- c('intercept', 'slope')

见 ?raster::calc