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

Question

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

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

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

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

我基于Function to sum each grid cells of raster stack using other rasters as an indicator 或此处Conditionally calculating the difference between max(raster) and each raster layer of raster stack 尝试了几种替代方案，但我不太了解它的逻辑。 我可以使用“字段计算器”使用 QGIS 解决同样的问题，并且它通过以下方式运行良好（它很乏味且可能出现错误）：

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

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

在哪里 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)