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

Question

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

我有 5 个 asciigrid 文件，代表 5 个光栅图像。有些像素确实有值，其他像素确实有 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 ... (at some point there will be numbers here)

编辑

感谢@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

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

Linear regression on raster images - lm complains about NAs

regression

r

raster

linear-regression

lm

编辑