CODEWITHSUNDEEP
rfunctiondplyrdata.table
89_Simple
89_Simple

Reputation: 3805

Applying a function to multiple groups in dplyr

I have a function that calculates soil water for a single day (from the ZeBook package)

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))
          }

this function takes three arguments: WAT0: water content of day i - 1 RAIN: Rainfall of day i,ETo: Evapotranspration of day i, CN, DC and FC which are constants.

It returns a dataframe with WAT1 which is the water content of day i, RO and DR

An example:

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

Now I want to run this function for day 1 till day 10. Sample data 

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

The below function uses the water.update function to calculate the soil water from day 1 till day 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)

This gives me three columns: first column with water content, second column is the RO and third is DR.

My issue is I need to run the 'water.model' function for multiple years and locations 

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))

I have two questions:

1) How do I run the water.model for each location and year from day 1 till day 10.

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

2) Any suggestions on making the above functions faster are welcome. Maybe using Rcpp? :)

EDIT

The function also takes a variable DC

Upvotes: 0

Views: 231

Answers (1)

Nate
Nate

Reputation: 10671

group_by() %>% nest() is a flexible tidyverse pattern for doing group-wise operations where the operation results can be different shapes/classes.

This example stores them as a list column, then you pull them out when you want with 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

Upvotes: 4

Related Questions