From e6a758d4fe19a3ba34d0f0dad1e373b5e0782f1c Mon Sep 17 00:00:00 2001 From: Beni Stocker Date: Tue, 29 Oct 2024 10:16:10 +0100 Subject: [PATCH 1/4] nothing really --- R/cwd.R | 2 -- R/simulate_snow.R | 2 -- 2 files changed, 4 deletions(-) diff --git a/R/cwd.R b/R/cwd.R index f30f92f..5d8b2f2 100644 --- a/R/cwd.R +++ b/R/cwd.R @@ -18,8 +18,6 @@ #' 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. #' -#' @importFrom dplyr -#' #' @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: diff --git a/R/simulate_snow.R b/R/simulate_snow.R index 28238a5..2b9a006 100644 --- a/R/simulate_snow.R +++ b/R/simulate_snow.R @@ -12,8 +12,6 @@ #' @param varnam_prec A character string specifying the variable name for rain. #' @param varnam_snow A character string specifying the variable name for snow. #' -#' @importFrom dplyr - #' @details Returns a data frame with two added columns: (1) \code{liquid_to_soil} #' is the rain plus snow melt in mm d-1; (2) \code{snow_pool} is the snow mass #' in water equivalents (mm) for each day. From 860e9759752114ee22da23673af3ddc6b1d9c885 Mon Sep 17 00:00:00 2001 From: Beni Stocker Date: Tue, 29 Oct 2024 10:53:01 +0100 Subject: [PATCH 2/4] implemented absolute threshold re-setting --- R/cwd.R | 16 ++++++++++++---- vignettes/potential_cwd.Rmd | 19 ++++++++++++++++--- 2 files changed, 28 insertions(+), 7 deletions(-) diff --git a/R/cwd.R b/R/cwd.R index 541cc71..bc6958d 100644 --- a/R/cwd.R +++ b/R/cwd.R @@ -13,6 +13,7 @@ #' 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. @@ -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 = 10, thresh_drop = 0.9, doy_reset = 999){ # corresponds to mct2.R @@ -69,12 +71,17 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d # continue accumulating deficit as long as the deficit has not fallen below (thresh_terminate) times the maximum deficit attained in this event # optionally while (iidx <= (nrow(df)-1) && # avoid going over row length - (deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit) + ((deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit) || + (deficit - df[[ varname_wbal ]][iidx] > max_deficit - thresh_terminate_absolute)) ){ dday <- dday + 1 deficit <- deficit - df[[ varname_wbal ]][iidx] + 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,13 +89,14 @@ 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_drop) && !done_finding_dropday || + 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 (df$doy[iidx] == doy_reset || deficit < (max_deficit - thresh_terminate_absolute)){ iidx_drop <- iidx max_deficit <- deficit break diff --git a/vignettes/potential_cwd.Rmd b/vignettes/potential_cwd.Rmd index a6144a0..36eaa99 100644 --- a/vignettes/potential_cwd.Rmd +++ b/vignettes/potential_cwd.Rmd @@ -182,6 +182,20 @@ 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)) ``` +Distributiono of daily precipitation. +```{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. ```{r} df <- df |> @@ -191,9 +205,8 @@ out_cwd <- cwd( df, varname_wbal = "pwbal", varname_date = "TIMESTAMP", - thresh_terminate = 0.0, - thresh_drop = 0.0, - doy_reset = doy_reset + thresh_terminate_absolute = 100.0, + thresh_drop = 0.0 ) ``` From 1ae5ff3a5251db444a833093a4b715e9a89362eb Mon Sep 17 00:00:00 2001 From: Patricia Helpap Date: Tue, 29 Oct 2024 16:47:12 +0100 Subject: [PATCH 3/4] changes implemented to include absolute cwd threshold and avoid negative values to be included in cwd events --- R/cwd.R | 13 ++++++++++--- vignettes/potential_cwd.Rmd | 12 +++++++++--- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/R/cwd.R b/R/cwd.R index bc6958d..0f2f25a 100644 --- a/R/cwd.R +++ b/R/cwd.R @@ -20,8 +20,8 @@ #' @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 @@ -69,8 +69,10 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, 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 >= 0) && # Ensure deficit is positive ############# changed to avoid negative deficit values ((deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit) || (deficit - df[[ varname_wbal ]][iidx] > max_deficit - thresh_terminate_absolute)) ){ @@ -78,6 +80,11 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, 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 + } + if (max_deficit > 0 && deficit < max_deficit - thresh_terminate_absolute){ print("now") } @@ -102,7 +109,7 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, 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 diff --git a/vignettes/potential_cwd.Rmd b/vignettes/potential_cwd.Rmd index 36eaa99..72700dd 100644 --- a/vignettes/potential_cwd.Rmd +++ b/vignettes/potential_cwd.Rmd @@ -182,7 +182,7 @@ 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)) ``` -Distributiono of daily precipitation. +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))) + @@ -197,6 +197,10 @@ df |> 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). + ```{r} df <- df |> mutate(pwbal = P_F - pet) @@ -205,6 +209,8 @@ out_cwd <- cwd( df, varname_wbal = "pwbal", varname_date = "TIMESTAMP", + #thresh_terminate = 2, + #doy_reset = 999, thresh_terminate_absolute = 100.0, thresh_drop = 0.0 ) @@ -221,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)), #reset to 0 for ylim 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)) + #reset to 0 for ylim labs( x = "Date", y = "Potential cumulative water deficit (mm)" From 93d4a28ed9f885052eaa8f751282b9680294500e Mon Sep 17 00:00:00 2001 From: Patricia Helpap Date: Wed, 30 Oct 2024 17:11:06 +0100 Subject: [PATCH 4/4] implemented changes to include an absolute termination threshold which defaults to NA --- R/cwd.R | 94 ++++++++++++++++++++++++++++++++----- vignettes/potential_cwd.Rmd | 14 +++--- 2 files changed, 89 insertions(+), 19 deletions(-) diff --git a/R/cwd.R b/R/cwd.R index 0f2f25a..4bb5d85 100644 --- a/R/cwd.R +++ b/R/cwd.R @@ -34,7 +34,7 @@ #' @export #' cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, - thresh_terminate_absolute = 10, thresh_drop = 0.9, doy_reset = 999){ + thresh_terminate_absolute = NA, thresh_drop = 0.9, doy_reset = 999){ # corresponds to mct2.R @@ -71,11 +71,80 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, # 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 >= 0) && # Ensure deficit is positive ############# changed to avoid negative deficit values - ((deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit) || - (deficit - df[[ varname_wbal ]][iidx] > max_deficit - thresh_terminate_absolute)) - ){ + 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] @@ -85,9 +154,10 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, break; ## Exit the loop if deficit is no longer positive } - if (max_deficit > 0 && deficit < max_deficit - thresh_terminate_absolute){ - print("now") - } + ## 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){ @@ -96,14 +166,13 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, } # record the day when deficit falls below (thresh_drop) times the current maximum deficit - if (deficit < (max_deficit * thresh_drop) && !done_finding_dropday || - deficit < (max_deficit - thresh_terminate_absolute) && !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 || deficit < (max_deficit - thresh_terminate_absolute)){ + if (deficit < (max_deficit - thresh_terminate_absolute)){ iidx_drop <- iidx max_deficit <- deficit break @@ -133,6 +202,7 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, df$deficit[iidx] <- deficit iidx <- iidx + 1 + } } diff --git a/vignettes/potential_cwd.Rmd b/vignettes/potential_cwd.Rmd index 72700dd..a28626c 100644 --- a/vignettes/potential_cwd.Rmd +++ b/vignettes/potential_cwd.Rmd @@ -199,7 +199,8 @@ df |> 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). +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 |> @@ -209,10 +210,9 @@ out_cwd <- cwd( df, varname_wbal = "pwbal", varname_date = "TIMESTAMP", - #thresh_terminate = 2, - #doy_reset = 999, - thresh_terminate_absolute = 100.0, - thresh_drop = 0.0 + #thresh_terminate = 0, + #thresh_terminate_absolute = 10.0, + #thresh_drop = 0.0 ) ``` @@ -227,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)), #reset to 0 for ylim + 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)) + #reset to 0 for ylim + ylim(0, max(out_cwd$df$deficit)) + labs( x = "Date", y = "Potential cumulative water deficit (mm)"