Merge pull request #8 from PHelpap/main

Include implementation of absolute threshold for cwd calculation

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,34 +69,116 @@ 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
           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){
+        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
+        }
 
       }
 

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