I am trying to work out a simple soil water balance in R. Here's the step I need to do:
For a given year, starting from doy 200,I need to determine the Soil Water (SWi) which is calculated by following formula
`SW(i) = SW(i-1) + Rain(i) - ETo(i)`
where SW(i-1) is the water content in previous day, Rain(i) is the rainfall and ETo(i) is the evapotranspiration on day i
The conditions are that SW(i) cannot be negative or be more than SW(max) which is 20.
Here's a sample data:
library(tidyverse)
set.seed(123)
dat <- tibble(
year = rep(1980:2015, each = 100),
day = rep(200:299, times = 36),
rain = sample(0:17, size = 100*36,replace = T),
eto = sample(2:9, size = 100*36,replace = T))
SW.initial <- data.frame(year= 1980:2015, SW.199 = sample(1:10, 36, replace = T))
SW.initial is the water content for doy 199 for for each year
SW.max <- 20
dat$SW.fin <- NA
Taking the example of year 1980
dat.1980 <- dat[dat$year == 1980,]
SW.initial.1980 <- SW.initial[SW.initial$year== 1980,"SW.199"]
for(doy in dat.1980$day){
SW <- SW.initial.1980
SW <- SW + dat.1980[dat.1980$day == doy, "rain"] - dat.1980[dat.1980$day == doy, "eto"]
SW <- ifelse(SW < 0, 0, ifelse(SW >= SW.max, SW.max, SW))
dat[dat$year == years & dat$day == doy,"SW.fin"] <- SW
SW.initial.1980 <- SW
}
This loop will give me the SW of each day starting doy 200 till 299 using:
`SW(i) = SW(i-1) + Rain[i] + ETo[i]``
where for doy 200, SW(i-1) was given from the SW.initial data frame
I can loop through all years:
for(years in unique(dat$year)){
test <- dat[dat$year == years,]
SW.in <- SW.initial[SW.initial$year == years,"SW.199"]
for(doy in test$day){
SW <- SW.in
SW <- SW + test[test$day == doy, "rain"] - test[test$day == doy, "eto"]
SW <- ifelse(SW < 0, 0, ifelse(SW >= SW.max, SW.max, SW))
dat[dat$year == years & dat$day == doy,"SW.fin"] <- SW
SW.in <- SW
}}
I really want to avoid this long loop and was thinking if there is much clever (and faster way to do this).
Does this give what you want ?
edit -> added grouped by year
dat %>% group_by(year) %>% mutate(sw_oneless = c(NA, day[1:length(day)-1]),
sw_oneless + rain - eto)
# A tibble: 3,600 x 6
# Groups: year [36]
year day rain eto sw_oneless `sw_oneless + rain - eto`
<int> <int> <int> <int> <int> <int>
1 1980 200 5 2 NA NA
2 1980 201 14 6 200 208
3 1980 202 7 4 201 204
4 1980 203 15 5 202 212
5 1980 204 16 5 203 214
6 1980 205 0 8 204 196
7 1980 206 9 9 205 205
8 1980 207 16 6 206 216
9 1980 208 9 9 207 207
10 1980 209 8 4 208 212
# ... with 3,590 more rows
To solve the "problem" with day 200, why don't you just filter out day 199-300 from your original data? You can then run my code and na.omit() or filter again and the rows with day 199 are gone. Or, if you can't do that, merge your SW.initial with your dat data frame