R中两个栅格图像之间的线性回归

Question

我需要进行线性回归，以计算一个经验参数。L1 是一种栅格图像，格式为 .tif。同时，L2 也是一种栅格图像，在先前已经计算过了。这两个图像具有相同的列数和行数。

公式如下： L1 = a + b * L2， 在 R 中表示为：

lm(L1 ~ L2)

在后面的公式中，我需要用到 a 和 b。

我现在遇到的问题是，这两个栅格都包含 NA 值，我不确定如何构建线性回归函数。我对 R 并不是那么熟悉，所以我被卡住了，而这个问题可能相当简单。我想我需要使用 calc，但不确定怎么做。 编辑：到目前为止，我有这段代码：

s = stack(L1,L2)
fun = function(x) {if (is.na(x[1])) { NA } else {lm(x[1] ~ x[2])$coefficients[2]}}

然而，计算时间非常长，且没有得出结果。

- Aldi Kasse 2

Answer 1

如果您想进行本地回归，即每个网格单元（像素）进行单独的回归，则应使用 calc。但在这种情况下这是没有意义的，因为您只有两个栅格；因此每个网格单元只有一个数据点。

在您的情况下，您似乎想要进行一次全局回归。您可以按照以下方法进行：

s <- stack(L1, L2)
v <- data.frame(na.omit(values(s)))
# this step may not be necessary
names(v) <- c('L1', 'L2')
m <- lm(L2 ~ L1, data=v)
m

如果s太大了，你可以尝试这样做

v <- sampleRegular(s, 100000)  
v <- data.frame(na.omit(v))

现在有了一些数据（并展示如何获取残差）

library(raster)
f <- stack(system.file("external/rlogo.grd", package="raster")) 
s <- stack(f)
names(s)
v <- data.frame(na.omit(values(s)))
m <- lm(red ~ green, data=v)
m

p <- predict(s, m)
residuals <- s$red - p

- Robert Hijmans

---

谢谢！我该如何添加一个条件，比如说：L2>0.0？如果这个条件不成立，就会出现逻辑错误的提示。 - Aldi Kasse 2

你可以对 v 进行子集操作。例如：v <- v[v$L2 > 0.0, ] - Robert Hijmans

你如何将这个模型的残差绘制成一个新的光栅图层？ - philiporlando

那将是一个新问题，但我已经在上面添加了它。 - Robert Hijmans