在dplyr中分组后运行函数

我编写了一个函数,用于计算作物生长阶段的时间(生长度天数)。作为背景,作物种植后,在积累了一定的热量单位后,会从一个阶段移动到另一个阶段。例如,对于给定的植物日,作物需要300°C、500°C、600°C的累积热量单位才能分别达到阶段1、阶段2和阶段3。

该函数取温度向量, temp.vec ,则, plant.date 这基本上是一天,你要开始计算累积热量单位,一个基本温度,一个最佳温度和一个临界温度。

set.seed(123)
sample.temp <- data.frame(day = 1:365,tmean = c(sample(25:32,365, replace = T)))

gdd.func <- function(temp.vec,plant.date,t.base,t.opt,t.cri){

 x <- temp.vec[temp.vec > plant.date]

 fT <- ifelse(x >= t.base & x <= t.opt,(x - t.base)/(t.opt - t.base),
           ifelse(t.opt <= x & x <= t.cri,(t.cri - x)/(t.cri - t.opt),0))

 Te <- t.base + fT*(t.opt - t.base)
 thermal.units <- Te - t.base
 day.stage1 <- which.max(cumsum(thermal.units) >= 300) # this will give me the day when cumulative accumulation of thermal units crossed 300 heat units

 # once growth stage 1 is reached, t.base,t.opt and t.cri are updated
 t.base <- t.base - 2
 t.opt <- t.opt - 2
 t.cri <- t.cri - 2

 fT[(day.stage1 + 1):length(fT)] <- ifelse(x[(day.stage1 + 1):length(fT)] >= t.base & x[(day.stage1 + 1):length(fT)] <= t.opt,(x[(day.stage1 + 1):length(fT)] - t.base)/(t.opt - t.base), ifelse(t.opt <= x[(day.stage1 + 1):length(fT)] & x[(day.stage1 + 1):length(fT)] <= t.cri,(t.cri - x[(day.stage1 + 1):length(fT)])/(t.cri - t.opt),0))

 Te[(day.stage1 + 1):length(Te)] <- t.base + fT[(day.stage1 + 1):length(fT)]*(t.opt - t.base)

 thermal.units[(day.stage1 + 1):length(Te)] <- Te[(day.stage1 + 1):length(Te)] - t.base

  day.stage2 <- which.max(cumsum(thermal.units) >= 500) 

   # once growth stage 2 is reached, t.base,t.opt and t.cri are updated again
   t.base <- t.base - 1
   t.opt <- t.opt - 1
   t.cri <- t.opt - 1

   fT[(day.stage2 + 1):length(fT)] <- ifelse(x[(day.stage2 + 1):length(fT)] >= t.base & x[(day.stage2 + 1):length(fT)] <= t.opt,(x[(day.stage2 + 1):length(fT)] - t.base)/(t.opt - t.base), ifelse(t.opt <= x[(day.stage2 + 1):length(fT)] & x[(day.stage2 + 1):length(fT)] <= t.cri,(t.cri - x[(day.stage2 + 1):length(fT)])/(t.cri - t.opt),0))

   Te[(day.stage2 + 1):length(Te)] <- t.base + fT[(day.stage2 + 1):length(fT)]*(t.opt - t.base)

   thermal.units[(day.stage2 + 1):length(Te)] <- Te[(day.stage2 + 1):length(Te)] - t.base

    day.stage3 <- which.max(cumsum(thermal.units) >= 600) 

    list(day.stage1,day.stage2,day.stage3)
   }

进行测试运行

   t.base <- 24
   t.opt <- 32
   t.cri <- 36

   plant.dates <- gdd.func(temp.vec = sample.temp$tmean,plant.date = 10,t.base,t.opt,t.cri)
   unlist(plant.dates)
  # [1]  66 117 144

输出是一个三天的向量,它给出了 植物日期 10

我的问题是,我是否要对多个 植物日期 跨多个地点和年份。例如,想象一下这些数据:

   sample.data <- data.frame(id1 = rep(1:20, each = 730*36), year = rep(rep(1980:2015, each = 365*2), times = 20),day = rep(rep(1:730, times = 36), times = 20), tmean = sample(25:32,20*730*36,replace = T))  

   head(sample.data)
   id1 year day tmean
1   1 1980   1    26
2   1 1980   2    32
3   1 1980   3    25
4   1 1980   4    26
5   1 1980   5    28
6   1 1980   6    28

数据由20个位置组成,每个位置有36年的数据。每年有730天(365*2)和每天的平均温度。

我有三个 植物日期 。

     plant.vec <- c(250,290,302)

我想选择每个种植日,并为我的每个位置X年组合生成三个生长阶段

    for(p in seq_along(plant.vec))

        plant.date <- plant.vec[p]

        sample.data %>% group_by(id1,year) %>% # how to insert my gdd.func here so that it runs for each id1 and year combination)

非常感谢。