Included re-setting CWD by day-of-year

DESCRIPTION

@@ -28,7 +28,9 @@ Suggests:
     ggplot2,
     patchwork,
     extRemes,
-    tidyr
+    tidyr,
+    visdat,
+    recipes
 License: AGPL-3
 RoxygenNote: 7.3.1
 Roxygen: list(markdown = TRUE)

R/cwd.R

@@ -34,7 +34,7 @@
 #'
 #' @export
 #'
-cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_drop = 0.9, doy_reset = NA){
+cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_drop = 0.9, doy_reset = 999){
 
   # corresponds to mct2.R
 
@@ -43,6 +43,9 @@ cwd <- function(df, varname_wbal, varname_date, thresh_terminate = 0.0, thresh_d
     thresh_terminate <- thresh_drop
   }
 
+  # create day-of-year column
+  df$doy <- lubridate::yday(df[[ varname_date ]])
+
   inst <- tibble()
   idx <- 0
   iinst <- 1
@@ -66,8 +69,10 @@ 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
-      while (iidx <= (nrow(df)-1) &&
-             (deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit)
+      # optionally
+      while (iidx <= (nrow(df)-1) &&  # avoid going over row length
+             (deficit - df[[ varname_wbal ]][iidx] > thresh_terminate * max_deficit) &&
+             df$doy[iidx] < doy_reset
              ){
 
         dday <- dday + 1

data-raw/get_data.R

@@ -13,6 +13,8 @@ saveRDS(df, file = paste0(here(), "/data/df_ch-lae.rds"))
 
 # get fluxnet data for FR-Pue site
 df <- read_csv(paste0(here(), "/inst/extdata/FLX_FR-Pue_FLUXDATAKIT_FULLSET_DD_2000_2014_2-3.csv")) |>
-  select(TIMESTAMP, P_F, TA_F_MDS, PA_F, LE_F_MDS, NETRAD)
+  select(TIMESTAMP, P_F, TA_F_MDS, PA_F, LE_F_MDS, SW_IN_F_MDS, LW_IN_F_MDS, NETRAD)
+
+# visdat::vis_miss(df)
 
 saveRDS(df, file = paste0(here(), "/data/df_fr-pue.rds"))

docs/pkgdown.yml

@@ -4,5 +4,5 @@ pkgdown_sha: ~
 articles:
   cwd_example: cwd_example.html
   potential_cwd: potential_cwd.html
-last_built: 2024-05-07T13:22Z
+last_built: 2024-05-13T13:11Z
 

vignettes/cwd_example.Rmd

@@ -110,11 +110,13 @@ gg5 / gg6 / gg7
 
 Get CWD and events.
 ```{r}
-out_cwd <- cwd(df,
-               varname_wbal = "wbal",
-               varname_date = "TIMESTAMP",
-               thresh_terminate = 0.0,
-               thresh_drop = 0.0)
+out_cwd <- cwd(
+  df,
+  varname_wbal = "wbal",
+  varname_date = "TIMESTAMP",
+  thresh_terminate = 0.0,
+  thresh_drop = 0.0
+  )
 ```
 
 Retain only events of a minimum length of 20 days.

vignettes/potential_cwd.Rmd

@@ -16,12 +16,15 @@ knitr::opts_chunk$set(echo = TRUE)
 ```{r message=FALSE}
 library(readr)
 library(dplyr)
+library(tidyr)
 library(here)
 library(lubridate)
 library(patchwork)
 library(extRemes)
 library(ggplot2)
 library(cwd)
+library(visdat)
+library(recipes)
 ```
 
 A potential cumulative water deficit can be calculated using net radiation and the *potential* evapotranspiration (PET). Here, we calculate PET based on net radiation and the approach by Priestley & Taylor (1972) as implemented by Davis et al. (2017). Necessary functions are part of the {cwd} package.
@@ -31,8 +34,60 @@ A potential cumulative water deficit can be calculated using net radiation and t
 Read data from the example FLUNXET file contained in this repository.
 ```{r}
 df <- readRDS(file = paste0(here(), "/data/df_fr-pue.rds"))
+
+visdat::vis_miss(df)
+```
+
+Some net radiation data is missing. Impute missing values by KNN.
+```{r}
+pp <- recipes::recipe(
+  NETRAD ~ SW_IN_F_MDS + LW_IN_F_MDS + TA_F_MDS,
+  data = df |> 
+    drop_na(SW_IN_F_MDS, LW_IN_F_MDS, TA_F_MDS)
+) |>
+  recipes::step_center(
+    recipes::all_numeric(),
+    -recipes::all_outcomes()
+  ) |>
+  recipes::step_scale(
+    recipes::all_numeric(),
+    -recipes::all_outcomes()
+  ) |>
+  recipes::step_impute_knn(
+    recipes::all_outcomes(),
+    neighbors = 5
+  )
+
+pp_prep <- recipes::prep(
+  pp,
+  training = df |> 
+    drop_na(SW_IN_F_MDS, LW_IN_F_MDS, TA_F_MDS)
+)
+
+df_baked <- recipes::bake(
+  pp_prep,
+  new_data = df
+)
+
+# fill missing with gap-filled
+df <- df |>
+  dplyr::bind_cols(
+    df_baked |>
+      dplyr::select(
+        NETRAD_filled = NETRAD)
+  ) |>
+  dplyr::mutate(
+    NETRAD = ifelse(is.na(NETRAD), NETRAD_filled, NETRAD)
+    #qc = ifelse(is.na(netrad), TRUE, FALSE)
+  ) |>
+  dplyr::select(
+    -NETRAD_filled
+  )
+
+visdat::vis_miss(df)
 ```
 
+
 ## Apply PET function
 
 ... and convert from units of mm s-1 to mm d-1.
@@ -68,8 +123,99 @@ df |>
   ) |> 
   ggplot() +
   geom_line(aes(doy, et, color = "ET")) +
-  geom_line(aes(doy, pet, color = "PET"))
+  geom_line(aes(doy, pet, color = "PET")) + 
+  labs(
+    x = "Day of year",
+    y = expression(paste("Water vapour mass flux (mm d"^-1, ")"))
+  ) +
+  theme_classic()
+```
+## Cumulating PET - *P*
+
+Check annual totals.
+```{r}
+adf <- df |>
+  mutate(year = year(TIMESTAMP)) |>
+  group_by(year) |>
+  summarise(pet = sum(pet), prec = sum(P_F))
+
+adf |>
+  tidyr::pivot_longer(cols = c(pet, prec), names_to = "Flux") |> 
+  ggplot(aes(x = year, y = value, color = Flux)) +
+  geom_line() +
+  labs(y = "Flux (mm/yr)") +
+  theme_classic()
+```
+
+In some cases, the mean annual PET may be larger than the mean annual precipitation (*P*), leading to a steady long-term increase of a *potential* cumulative water deficit. This is not the case here (see plot above). For demonstration, let's assume precipitation was 30% of its actual value and calculate the running sum of (PET - *P*) - the cumulative potential evapotranspiration.
+```{r}
+df |> 
+  mutate(P_F = 0.3 * P_F) |>
+  mutate(pcwd = cumsum(pet - P_F)) |> 
+  ggplot(aes(TIMESTAMP, pcwd)) +
+  geom_line() +
+  theme_classic()
+```
 
+This indicates a need to re-set the potential cumulative water deficit calculation. Let's determine the wettest month from the available years of data and reset the cumulative water deficit each year in that month. The plot below shows the average *P* - PET for each month. November is the wettest month at this site. 
+
+```{r}
+mdf_mean <- df |> 
+  mutate(month = lubridate::month(TIMESTAMP),
+         pwbal = P_F - pet) |> 
+  group_by(month) |> 
+  summarise(
+    pwbal = sum(pwbal)
+  )
+
+mdf_mean |> 
+  ggplot(aes(as.factor(month), pwbal)) +
+  geom_bar(stat = "identity") +
+  theme_classic()
+```
+
+Therefore, we re-set the cumulation of the water deficit each year in November.
+```{r}
+# determine the day-of-year of the first day of the month after the wettest month
+doy_reset <- lubridate::yday(lubridate::ymd("2000-11-01") + lubridate::dmonths(1))
+```
+
+Get potential CWD and events.
+```{r}
+df <- df |> 
+  mutate(pwbal = P_F - pet)
+
+out_cwd <- cwd(
+  df,
+  varname_wbal = "pwbal",
+  varname_date = "TIMESTAMP",
+  thresh_terminate = 0.0,
+  thresh_drop = 0.0,
+  doy_reset = doy_reset
+  )
+```
+
+Retain only events of a minimum length of 20 days.
+```{r}
+out_cwd$inst <- out_cwd$inst |>
+  filter(len >= 20)
+```
+
+Plot CWD time series.
+```{r}
+ggplot() +
+  geom_rect(
+    data = out_cwd$inst,
+    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)) +
+  labs(
+    x = "Date", 
+    y = "Potential cumulative water deficit (mm)"
+    )
 ```
 
 ## References