如何基于Raster Brick/Stack创建列值超阈值次数的栅格图层

现有如下R代码示例，创建了一个简单的Raster Brick，并通过calc函数生成了记录每个像元最大值所在图层编号的栅格图层。现在需要修改代码，生成一个显示每个像元数值超过指定值（例如1）的次数的栅格图层。

library(raster)

## generate simple raster brick
r1 <- r2 <- r3 <- r4 <- r5 <- raster(ncol=10, nrow=10)
r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
r4[] <- rnorm(ncell(r4), 1, 0.2)
r5[] <- rnorm(ncell(r5), 1, 0.2)

s <- brick(r1,r2,r3,r4,r5)

## find which layer the max value occurs
s_max <- calc(s, which.max)

plot(s_max)

解决方案 #

只需在calc函数中传入一个自定义匿名函数，统计每个像元跨图层数值中大于指定阈值的元素数量即可。修改后的代码如下：

library(raster)

## generate simple raster brick
r1 <- r2 <- r3 <- r4 <- r5 <- raster(ncol=10, nrow=10)
r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
r4[] <- rnorm(ncell(r4), 1, 0.2)
r5[] <- rnorm(ncell(r5), 1, 0.2)

s <- brick(r1,r2,r3,r4,r5)

## 统计每个像元数值超过1的次数
count_over_threshold <- calc(s, function(x) sum(x > 1, na.rm = TRUE))

plot(count_over_threshold)

代码说明 #

• function(x) sum(x > 1, na.rm = TRUE)：对每个像元的跨图层数值（x是包含该像元所有图层值的向量），先判断每个值是否大于1（布尔值会自动转换为1/0），再求和得到次数。
• na.rm = TRUE：自动忽略栅格中的NA值，避免统计结果出现NA。若无需处理NA，可移除该参数。
• 调整阈值：直接将x > 1中的1替换为目标数值即可。

内容的提问来源于stack exchange，提问作者Blaiso