通过 R 中的栅格堆栈条件求和

Question

我有一个 34 组光栅 (.tif)，每个光栅大约 10MB，值 1 和 0 分别代表桉树的存在或不存在，相同区域和相同的监督分类的产物。每个栅格代表从1985年到2018年的一年评估（i）。我想在三个图像组中找到两种类型的像素序列的总和：

稳定像素：第(i)年像素值为1，(i+1) e(i-1)对应一个111的序列，也就是说连续多年是桉树。
不稳定像素：1英寸第 i 年和第 i + 1 年和第 i-1 年（序列 010）中的 0。那就是说它交替是桉树，而在其他年份则不是，这可能是分类错误。

在这两种情况下，目的都是验证分类是否稳健，是否超出分类产生的高 kappa 指数。我的问题是：

如何在 R 中实现这个过程？
我读到 "for" 的循环可能不是处理大型栅格的最佳选择，还有其他选择吗？因为我要实现类似逻辑的其他流程。

我尝试了几种基于： or here 的替代方案，但我不太了解它的逻辑。我可以通过 "Field Calculator" 使用 QGIS 解决同样的问题，并且它通过以下方式运行良好（它很乏味且容易出错）：

稳定像素 = ("img1985 @ 1" = 1 AND "img1986 @ 1" = 1 AND "img1987 @ 1" = 1) + ... + ("img2016 @ 1" = 1 AND "img2017 @ 1" = 1 和 "img2018 @ 1" = 1)

不稳定像素 = ("img1985 @ 1" = 0 AND "img1986 @ 1" = 1 AND "img1987 @ 1" = 0) + ... + ("img2016 @ 1" = 0 AND "img2017 @ 1" = 1 和 "img2018 @ 1" = 0)

where to imgYear: 是每一年的每个栅格； 1985 ... 2018: 年

```
library(raster)
# Initial sample data + result
r1 <- r2 <- r3 <- r4 <- r5 <- rUns <- rSta <- raster(matrix(0, 10, 10))

# Create some "stable pixels" of example.
for (i in seq(from=10, to=50, by=10)) {
r1[(i+6):(i+10)] <- 1
r2[(i+6):(i+10)] <- 1
r3[(i+6):(i+10)] <- 1
r4[(i+6):(i+10)] <- 1
r5[(i+6):(i+10)] <- 1
}
# Create some "unstable pixels" of example.
r2 [c(60,70,80)] <- 1
r3 [1] <- 1
r4 [c(60,70,80)] <- 1
# Stack raster
r <- stack(r1, r2, r3, r4, r5)

# *** Expected results ***
# Sum of stable pixels (Sta)
for (i in 2:nlayers(r)) {
pixelSta <- ((r[[i-1]] == 1) * (r[[i]] == 1) * (r[[i+1]] == 1))* 1
}
# Don't work: Error in .local(x, ...) : not a valid subset
# *** Result Expected (rSta): sequence 111. Manually.
for (i in seq(from=10, to=50, by=10)) {
rSta[(i+6):(i+10)] <- 3
}
as.matrix(rSta)

# Sum of unstable pixels (Uns)
for (i in 2:nlayers(r)) {
pixelUns <- ((r[[i-1]] == 0) * (r[[i]] == 1) * (r[[i+1]] == 0))* 1
}
# Don't work: Error in .local(x, ...) : not a valid subset
# *** Result Expected (rUns): sequence 010. Manually.
rUns [c(60,70,80)] <- 2
rUns [1] <- 1
as.matrix(rUns)
```

预期结果在代码中（*** 预期结果，手动）。`

非常感谢你，我希望我已经说清楚了。

Answer 1

您需要一个计算稳定性的简单函数。例如，函数 f

f <- function(x) sum(diff(x)==0)

计算不是从一种状态转换到另一种状态的步数。

试一试

f(c(0,0,0,0,0,0))
#[1] 5
f(c(0,0,0,0,1,1))
#[1] 4
f(c(0,1,0,1,0,1))
#[1] 0

现在使用栅格数据

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster"))
s <- stack(s, s)
x <- calc(s, f)
plot(x)

通过 R 中的栅格堆栈条件求和

Sum conditioned through raster stack in R

r

raster

r-raster