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

问题描述

我确定可以用几个字节来解决这个问题，但是我已经花了数小时来处理这个简单的事情，并且无法摆脱.我不经常使用R.

I'm sure this can be fixed with few bytes, but I've spent hours on this simple thing and can't get out of it. I don't use R often.

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

I have 5 asciigrid files that represent 5 raster images. Some pixels do have values, other do have NAs. For example, the first image might be something like:

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％确信.我愿意做的事:

As you can see, NA position is always the same and I'm 100% sure about that. What I'm willing to do:

• 使用 read.asciigrid();
读取文件
• 使用 raster 包中的 values()在长数组中获取它们的值；
• 创建一个包含5行的矩阵，每行包含相应地图的值；
• 线性拟合各列并获得系数.每列将代表一个像素，并具有与5张地图相对应的5个值.
• 使用系数值创建两个新的光栅图像.

• read the files with read.asciigrid();
• get their values in a long array, using values() from the raster package;
• create a matrix with 5 rows, each of them holding the values of the corresponding map;
• linear fit each column and get the coefficients. Each column will represent a pixel, and will have 5 values corresponding to the 5 maps.
• create two new raster images with the coefficient values.

我被困在 lm .具体来说，它说: lm.fit(...)中的错误:0(非NA)情况.但是，据我了解的归因贴图，应该要么具有 all NA的列，要么完全具有 no NA的列，如下所示:

I'm stuck at lm. Specifically, it says: Error in lm.fit(...): 0 (non-NA) cases. However, from what I know about the imput maps, there should be either columns with all NAs or columns with no NA at all, like as follows:

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位置.我哪里错了?在下面粘贴我的代码.谢谢.

so I can create a new raster image with the coefficients and keep the NA position. Where am I wrong? Pasting my code below. Thank you.

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)的输出:

编辑

由于@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 调用上使用此功能.我看着 raster :::.calcTest 没有成功.我尝试管理所有 y 值均为 NA 的情况，如下所示:

Getting an obscure Cannot use this function on a .calcTest call. I looked at raster:::.calcTest with no success. I have tried managing the case where all y values are NA, like so:

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

现在它工作了几分钟，但后来我得到 Error in setValues(out, x) :值必须是数字、整数、逻辑或因子.但是，通常将NA设置为栅格值！我看不到这是怎么回事.

And now It works for some minutes, but then I'm getting Error in setValues(out, x) : values must be numeric, integer, logical or factor. However it is common to set NA to raster values! I can't see what's wrong here.

推荐答案

这是获取栅格数据的方法

This is how you can get the raster data

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

您现在可以执行类似的操作.

You could now do something like this.

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

但是您不需要这样做.为了做这样的回归，我会做

But you do not need to do that. To do regression like this, I would do

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

但是，因为您的单元格(各层之间)仅具有NA值，所以这将失败(如上述应用中所述).您需要编写一个函数来捕获这些情况:

But because you have cells with only NA values (across the layers) this will fail (like in the apply above). You need to write a function to catch these cases:

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

See ?raster::calc

这篇关于光栅图像上的线性回归-LM抱怨NA的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！