如何计算 R 中使用两个栅格堆栈的每个组的最大值？答案

【问题标题】：How to calculate Max value of each group used two raster stack in R?如何计算 R 中使用两个栅格堆栈的每个组的最大值？
【发布时间】：2021-01-15 01:50:18
【问题描述】：

我想根据另一个光栅堆栈的组对光栅堆栈（包含多个图层）应用计算（总和）（组标识从一个像素时间序列到另一个像素时间序列是可变的）。这是我的例子：

library(raster)
#1. Create raster stack "values" (rs)
set.seed(193473)
r1 <- raster(nrows = 10, ncols = 10)
r2=r3=r4=r5=r6=r7=r8=r9=r10=r11=r12=r13=r1

r1[] <- rbinom(ncell(r1), 60, prob = .1)
r2[] <- rbinom(ncell(r1), 50, prob = .1)
r3[] <- rbinom(ncell(r1), 70, prob = .1)
r4[] <- rbinom(ncell(r1), 90, prob = .1)

r5[] <- rbinom(ncell(r1), 70, prob = .1)
r6[] <- rbinom(ncell(r1), 60, prob = .1)
r7[] <- rbinom(ncell(r1), 80, prob = .1)
r8[] <- rbinom(ncell(r1), 90, prob = .1)

r9[] <- rbinom(ncell(r1), 60, prob = .1)
r10[] <- rbinom(ncell(r1), 70, prob = .1)
r11[] <- rbinom(ncell(r1), 90, prob = .1)
r12[] <- rbinom(ncell(r1), 70, prob = .1)
r13[] <- rbinom(ncell(r1), 50, prob = .1)
  
rs <- stack(r1, r2, r3, r4,r5,r6,r7,r8,r9,r10,r11,r12,r13)
nlayers(rs)
plot(rs[[1]])
plot(rs)

#2, rs_flag

rc1=function(x1) {
  ifelse(x1>7,1,0)
}

rs_flag=overlay(rs,fun=rc1)
nlayers(rs_flag)
plot(rs_flag)

#3,  rs_id

rc3= function(x3) {
  rep(seq(1,length(rle(x3)$lengths)), rle(x3)$lengths)
}

rs_id=overlay(rs_flag,fun=rc3)

plot(rs_id)

如何提取每组光栅栈(rs_id)中(rs)光栅栈的最大值，并用这个最大值替换同组的其他值？

我解释了我的问题，请看下面一个像素的时间序列示例：

#extract values of one pixel ex. (x,y)=(4,7)
rsp=rs[4,7]
rsp

#rsp_id
rsp_id=rs_id[4,7]
rsp_id

提取每组堆栈栅格（rs_id）的最大值，并将同一组的其他值替换为该最大值。

#table preparation  
tab1 <- t(rbind(rsp, rsp_id))
tab1

#rsp_max : max value of each group

library(dplyr)
   
tab2=tab1 %>%
     group_by(rsp_id) %>%
     mutate(rsp_max=max(rsp))

tab2

最后我得到一个像素的结果（见下表）

我的问题：如何在函数内部使用“group_by”和“max”在两个光栅堆栈（rs 和 rs_id）之间应用上述函数来计算“rs_max”（光栅堆栈）

#我试试这个功能但是不行！！：

##4# rs_id

library(dplyr)

rc4= function(x,y) {
  group_by(x) %>%
   max(y)
}



#Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE,  : 
#cannot use this formula, probably because it is not vectorized

提前感谢您的帮助！

【问题讨论】：

标签： r dplyr statistics r-raster

【解决方案1】：

我认为这可以解决问题

ff <- function(v) {
    x <- v[1:13]
    y <- v[14:26]
    a <- tapply(x, y, max) 
    b <- as.integer(names(a))
    as.vector(a[match(y, b)])
}

ff(as.vector(tab1))

rs_max <- calc(stack(rs, rs_id), fun=ff)
rs_max[4,7]

【讨论】：

非常感谢 Robert Hijmans 的回答。该功能运行良好。感谢这个想法：堆叠两个栅格堆栈并使用tapply。