【发布时间】:2018-03-27 00:16:45
【问题描述】:
我编写了一个函数来计算作物生长阶段的时间(生长期天数)。作为背景,作物种植后,在积累了一定的热量单位后,从一个阶段转移到另一个阶段。例如对于给定的种植日,作物需要 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
输出是一个三天的向量,它给出了plant.date 10 的 stage1、stage2 和 stage3。
我的问题是我是否想跨多个位置和年份为多个plant.date 运行上述函数。例如想象一下这个数据:
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.date。
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)
谢谢
【问题讨论】:
-
我建议如果你大大减少这个问题,你会得到更多的回应。它相当长,而且很可能 80-90% 的代码不是问题的一个因素。我个人觉得这个长度有点太吓人了,无法快速浏览一下,一旦它从我的视线中消失......
-
好的..让我现在缩短这个。
标签: r dplyr data.table apply purrr