Применение функции к нескольким группам в 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) Как запустить water.model для каждого местоположения и года с 1 по 10 день.

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