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

编程问答 2022-06-11

问题描述

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

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),prob = .1)
r6[] <- rbinom(ncell(r1),prob = .1)
r7[] <- rbinom(ncell(r1),80,prob = .1)
r8[] <- rbinom(ncell(r1),prob = .1)

r9[] <- rbinom(ncell(r1),prob = .1)
r10[] <- rbinom(ncell(r1),prob = .1)
r11[] <- rbinom(ncell(r1),prob = .1)
r12[] <- rbinom(ncell(r1),prob = .1)
r13[] <- rbinom(ncell(r1),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

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

enter image description here

我的问题：如何在函数内部使用“ 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

提前感谢您的帮助！

解决方法

我认为这可以解决问题

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]
dplyrdplyrrrr-rasterstatistics