daily pastebin goal
63%
SHARE
TWEET

Untitled

a guest May 16th, 2018 102 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. water.update <- function(WAT0, RAIN, ETo){
  2.  
  3.             S = 25400/CN - 254;  IA = 0.2*S
  4.  
  5.             if(RAIN > IA){RO = (RAIN - 0.2 * S)^2/(RAIN + 0.8 * S)
  6.             } else {
  7.             RO = 0
  8.             }
  9.  
  10.             if(WAT0 + RAIN - RO > FC) {DR = DC * (WAT0 + RAIN - RO - FC)
  11.             } else {
  12.             DR = 0
  13.             }    
  14.             dWAT = RAIN - RO - DR - ETo
  15.             WAT1 = WAT0 + dWAT
  16.             return(c(WAT1,RO,DR))
  17.           }
  18.    
  19. CN <- 60;FC <- 42;DC <- 0.02
  20. water.update(WAT0 = 23, RAIN = 5, ETo = 2)
  21. # 26, 0, 0
  22.    
  23. weather <- data.frame(day = 1:10 ,rain = sample(1:100, 10, replace = T), ETo = sample(1:10, 10, replace = T))
  24.    
  25. water.model <- function(weather, FC, CN, WAT0){
  26.  
  27.     WAT <- data.frame(matrix(NA, nrow =  nrow(weather), ncol = 3))
  28.     WAT[1,1] <- WAT0 # WAT0 is a constant
  29.  
  30.     for(day in 1:(nrow(weather)-1)){
  31.  
  32.       WAT[day + 1,] = water.update(WAT[day,1],weather$rain[day],weather$ETo[day])
  33.     }
  34.     return(WAT)
  35.     }
  36.  
  37. WAT0 <- 20
  38. water.model(weather = weather, FC = FC, CN = CN, WAT0 = WAT0)
  39.    
  40. big.data <- data.frame(loc.id = rep(1:3, each = 10*3),
  41.                             year = rep(rep(1981:1983, each = 10),times = 3),
  42.                             day = rep(1:10, times = 3*3),
  43.                             CN = rep(c(50,55,58), each = 10*3), # each location has a contant CN, FC and DC
  44.                             FC = rep(c(72,76,80),each = 10*3),
  45.                             DC = rep(c(0.02,0.5,0.8), each = 10*3),
  46.                             WAT0 = rep(c(20,22,26), each = 10*3),
  47.                             rain = sample(1:100,90, replace = T),
  48.                             eto = sample(1:10,90, replace = T))
  49.    
  50. big.data %>% group_by(loc.id, year) %>% do??
RAW Paste Data
Top