diff --git a/R/cwd.R b/R/cwd.R index 541cc71..4bb5d85 100644 --- a/R/cwd.R +++ b/R/cwd.R @@ -13,14 +13,15 @@ #' during the same event, to which a CWD has to be reduced to terminate the event. Defaults to 0, #' meaning that the CWD has to be fully compensated by water infiltration into the soil to terminate #' a CWD event. +#' @param thresh_terminate_absolute threshold determining end of event, as \code{thresh_terminate} but in absolute terms (in mm d-1). #' @param thresh_drop Level, relative to the CWD maximum of the same event, after which all data #' during the remainder of the event is set to missing values. This is to avoid interpreting data #' after rain events but before full compensation of CWD. Defaults to 0.9. #' @param doy_reset Day-of-year (integer) when deficit is to be reset to zero each year. #' #' @details A list of two data frames (tibbles). \code{inst} contains information about CWD "events". -#' Each row corresonds to one event. An event is defined as a period of consecutive days where the -#' CWD is positive (a water deficit)) and has the following columns: +#' Each row corresponds to one event. An event is defined as a period of consecutive days where the +#' CWD is positive (a water deficit) and has the following columns: #' #' \code{idx_start}: row number of \code{df} of which the date corresponds to the start of the event #' \code{len}: length of the event, quantified as number of rows in \code{df} corresponding to the event @@ -32,7 +33,8 @@ #' #' @export #' -cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_drop = 0.9, doy_reset = 999){ +cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, + thresh_terminate_absolute = NA, thresh_drop = 0.9, doy_reset = 999){ # corresponds to mct2.R @@ -67,14 +69,96 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d done_finding_dropday <- FALSE # continue accumulating deficit as long as the deficit has not fallen below (thresh_terminate) times the maximum deficit attained in this event + # OR as long as deficit has not fallen below the maximum deficit attained - (thresh_terminat_absolut) # optionally - while (iidx <= (nrow(df)-1) && # avoid going over row length - (deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit) - ){ + if(is.na(thresh_terminate_absolute)) {######determines whether an argument is given for thresh_terminate_absolute + + while (iidx <= (nrow(df)-1) && # avoid going over row length + (deficit >= 0) && # Ensure deficit is positive + ((deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit)) + ){ + + dday <- dday + 1 + deficit <- deficit - df[[ varname_wbal ]][iidx] + + ##Immediately stop if deficit falls below zero + if (deficit < 0) { + break; ## Exit the loop if deficit is no longer positive + } + + ## for development: + # if (max_deficit > 0 && deficit < max_deficit - thresh_terminate_absolute){ + # print("now") + # } + + # record the maximum deficit attained in this event + if (deficit > max_deficit){ + max_deficit <- deficit + done_finding_dropday <- FALSE + } + + # record the day when deficit falls below (thresh_drop) times the current maximum deficit + if (deficit < (max_deficit * thresh_drop) && !done_finding_dropday) { + iidx_drop <- iidx + done_finding_dropday <- TRUE + } + + # stop accumulating on re-set day + if (df$doy[iidx] == doy_reset) { + iidx_drop <- iidx + max_deficit <- deficit + break + } + + # once deficit has fallen below threshold, all subsequent dates are dropped (dday set to NA) + if (done_finding_dropday){ + df$iinst[iidx] <- NA + df$dday[iidx] <- NA + } else { + df$iinst[iidx] <- iinst + df$dday[iidx] <- dday + iidx_drop <- iidx + } + + # # (even before "drop-day" is found), drop data of days after rain, i.e., where current CWD is below the maximum (previously) attained in the same event + # if (deficit < max_deficit){ + # df$iinst[iidx] <- NA + # df$dday[iidx] <- NA + # df$deficit[iidx] <- NA + # } else { + # df$iinst[iidx] <- iinst + # df$dday[iidx] <- dday + # df$deficit[iidx] <- deficit + # } + + df$deficit[iidx] <- deficit + + iidx <- iidx + 1 + + + + } + + } else { ###########starts counting instances if thresh_terminate_absolute is given + + while (iidx <= (nrow(df)-1) && # avoid going over row length + (deficit >= 0) && # Ensure deficit is positive + ((deficit - df[[ varname_wbal ]][iidx] > max_deficit - thresh_terminate_absolute)) + ){ dday <- dday + 1 deficit <- deficit - df[[ varname_wbal ]][iidx] + ##Immediately stop if deficit falls below zero + if (deficit < 0) { + break; ## Exit the loop if deficit is no longer positive + } + + ## for development: + # if (max_deficit > 0 && deficit < max_deficit - thresh_terminate_absolute){ + # print("now") + # } + # record the maximum deficit attained in this event if (deficit > max_deficit){ max_deficit <- deficit @@ -82,19 +166,19 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d } # record the day when deficit falls below (thresh_drop) times the current maximum deficit - if (deficit < (max_deficit * thresh_drop) && !done_finding_dropday){ + if (deficit < (max_deficit - thresh_terminate_absolute) && !done_finding_dropday){ iidx_drop <- iidx done_finding_dropday <- TRUE } # stop accumulating on re-set day - if (df$doy[iidx] == doy_reset){ + if (deficit < (max_deficit - thresh_terminate_absolute)){ iidx_drop <- iidx max_deficit <- deficit break } - # once, deficit has fallen below threshold, all subsequent dates are dropped (dday set to NA) + # once deficit has fallen below threshold, all subsequent dates are dropped (dday set to NA) if (done_finding_dropday){ df$iinst[iidx] <- NA df$dday[iidx] <- NA @@ -118,6 +202,7 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d df$deficit[iidx] <- deficit iidx <- iidx + 1 + } } diff --git a/vignettes/potential_cwd.Rmd b/vignettes/potential_cwd.Rmd index a6144a0..a28626c 100644 --- a/vignettes/potential_cwd.Rmd +++ b/vignettes/potential_cwd.Rmd @@ -182,7 +182,26 @@ Therefore, we re-set the accumulation of the water deficit each year on the firs doy_reset <- lubridate::yday(lubridate::ymd("2000-11-01") + lubridate::dmonths(1)) ``` +Distribution of daily precipitation and PET to aid choosing an absolute threshold (`threshold_terminate_absolute`) determining the end of a CWD event (see next). +```{r} +df |> + ggplot(aes(P_F, after_stat(density))) + + geom_density(bw = 1) + + xlim(-1, 30) + +df |> + ggplot(aes(pet, after_stat(density))) + + geom_density(bw = 1) + + xlim(-1, 30) +``` + + Get the potential cumulative water deficit time series and individual *events*. Note that we use the argument `doy_reset` here to force a re-setting of the potential cumulative water deficit on that same day each year. +A relative and an absolute option are possible for determining the threshold to which a CWD has to be reduced to to terminate the event: +1) `threshold_terminate`: Set the threshold relative to maximum CWD attained. Defaults to 0, meaning that the CWD has to be fully compensated by water infiltration into the soil to terminate a CWD event. +2) `threshold_terminate_absolute`: threshold determining end of event as `thresh_terminate` but in absolute terms (in mm d-1). Default set to `NA`. +If no option is given, the function defaults to the relative threshold. + ```{r} df <- df |> mutate(pwbal = P_F - pet) @@ -191,9 +210,9 @@ out_cwd <- cwd( df, varname_wbal = "pwbal", varname_date = "TIMESTAMP", - thresh_terminate = 0.0, - thresh_drop = 0.0, - doy_reset = doy_reset + #thresh_terminate = 0, + #thresh_terminate_absolute = 10.0, + #thresh_drop = 0.0 ) ``` @@ -208,12 +227,12 @@ Plot the potential cumulative water deficit time series and events. ggplot() + geom_rect( data = out_cwd$inst, - aes(xmin = date_start, xmax = date_end, ymin = 0, ymax = max( out_cwd$df$deficit)), + aes(xmin = date_start, xmax = date_end, ymin = 0, ymax = max( out_cwd$df$deficit)), fill = rgb(0,0,0,0.3), color = NA) + geom_line(data = out_cwd$df, aes(TIMESTAMP, deficit), color = "tomato") + theme_classic() + - ylim(0, max( out_cwd$df$deficit)) + + ylim(0, max(out_cwd$df$deficit)) + labs( x = "Date", y = "Potential cumulative water deficit (mm)"