将函数应用于 dplyr 中的多个组

Question

我有一个计算一天土壤水分的函数（来自 ZeBook 包）

water.update <- function(WAT0, RAIN, ETo){

            S = 25400/CN - 254;  IA = 0.2*S

            if(RAIN > IA){RO = (RAIN - 0.2 * S)^2/(RAIN + 0.8 * S)
            } else { 
            RO = 0 
            }

            if(WAT0 + RAIN - RO > FC) {DR = DC * (WAT0 + RAIN - RO - FC) 
            } else { 
            DR = 0 
            }    
            dWAT = RAIN - RO - DR - ETo
            WAT1 = WAT0 + dWAT
            return(c(WAT1,RO,DR))
          } 

这个函数接受三个参数：WAT0: 第 i - 1 天的含水量 RAIN：第 i 天降雨，ETo：第 i 天蒸发量，CN，DC 和 FC 是常数。

它返回一个带有 WAT1 的数据帧，它是第 i 天、RO 和 DR 的含水量

一个例子：

CN <- 60;FC <- 42;DC <- 0.02
water.update(WAT0 = 23, RAIN = 5, ETo = 2)
# 26, 0, 0 

现在我想在第 1 天到第 10 天运行这个函数。示例数据

weather <- data.frame(day = 1:10 ,rain = sample(1:100, 10, replace = T), ETo = sample(1:10, 10, replace = T))

下面的函数使用water.update函数计算第1天到第10天的土壤水分。

water.model <- function(weather, FC, DC,CN, WAT0){

    WAT <- data.frame(matrix(NA, nrow =  nrow(weather), ncol = 3))
    WAT[1,1] <- WAT0 # WAT0 is a constant 

    for(day in 1:(nrow(weather)-1)){

      WAT[day + 1,] = water.update(WAT[day,1],weather$rain[day],weather$ETo[day])
    }
    return(WAT)
    }

WAT0 <- 20
water.model(weather = weather, FC = FC, CN = CN, WAT0 = WAT0)

这给了我三列：第一列是含水量，第二列是 RO，第三列是 DR。

我的问题是我需要在多年和多个地点运行“water.model”功能

big.data <- data.frame(loc.id = rep(1:3, each = 10*3), 
                            year = rep(rep(1981:1983, each = 10),times = 3), 
                            day = rep(1:10, times = 3*3),
                            CN = rep(c(50,55,58), each = 10*3), # each location has a contant CN, FC and DC
                            FC = rep(c(72,76,80),each = 10*3),
                            DC = rep(c(0.02,0.5,0.8), each = 10*3),
                            WAT0 = rep(c(20,22,26), each = 10*3),
                            rain = sample(1:100,90, replace = T),
                            eto = sample(1:10,90, replace = T))

我有两个问题：

1) 如何从第 1 天到第 10 天为每个位置和年份运行 water.model。

big.data %>% group_by(loc.id, year) %>% do??

2) 欢迎任何关于使上述功能更快的建议。也许使用 Rcpp？ :)

编辑

该函数还接受一个变量DC

Answer 1

group_by() %>% nest() 是一种灵活的tidyverse 模式，用于进行分组操作，其中操作结果可以是不同的形状/类。

此示例将它们存储为列表列，然后您可以在需要时使用unnest() 将它们拉出。

# redefine for typo in var name, your fxn expects 'ETo' not 'eto'
big.data <- data.frame(loc.id = rep(1:3, each = 10*3), 
                       year = rep(rep(1981:1983, each = 10),times = 3), 
                       day = rep(1:10, times = 3*3),
                       CN = rep(c(50,55,58), each = 10*3),
                       FC = rep(c(72,76,80),each = 10*3),
                       DC = rep(c(0.02,0.5,0.8), each = 10*3),
                       WAT0 = rep(c(20,22,26), each = 10*3),
                       rain = sample(1:100,90, replace = T),
                       ETo = sample(1:10,90, replace = T)) # typo was here

library(tidyverse)

res <- big.data %>%
  group_by(loc.id, year) %>%
  nest() %>%
  mutate(mod = map(data, ~ water.model(weather = ., FC = FC, CN = CN, WAT0 = WAT0)))

unnest(res, mod)

# A tibble: 90 x 5
   loc.id  year    X1    X2     X3
    <int> <int> <dbl> <dbl>  <dbl>
 1      1  1981  20.0 NA    NA    
 2      1  1981  78.5  6.68  0.846
 3      1  1981 124.   2.14  1.77 
 4      1  1981 190.  14.4   3.16 
 5      1  1981 238.   2.34  4.01 
 6      1  1981 297.  11.1   5.35 
 7      1  1981 357.  11.1   6.54 
 8      1  1981 349.   0.    6.37 
 9      1  1981 400.   6.35  7.42 
10      1  1981 420.   0.    7.81 
# ... with 80 more rows