R/fExtract_SOSEOS_bySite.R

In ksmiff33/FluxSynthU: Suite of functions to analyze eddy covariance data

#' This function pulls out critical threshold dates (SOS/EOS) from 466-day time series of GPPsat.

#' @export
#' @title Estimate critical dates (SOS/EOS) from GPPsat time series.

#' @param POSIX date, vector of POSIX dates (format YYYY-MM-DD)
#' @param cumu.DoY numeric, vector of cumulative day of year based on a grouping interval
#' @param proc_yr numeric, vector of process years by which to group analysis
#' @param GPP numeric, vector of GPP data for determining SOS/EOS (e.g. GPPsat)



fExtract_SOSEOS_bySite <- function(POSIX, cumu.DoY, proc_yr, GPP) {

  POSIX <- as.Date(POSIX, "%Y-%m-%d")
  DoY <- yday(POSIX)

  df <- data.frame(POSIX, cumu.DoY, proc_yr, GPP)
  names(df) <- c("k_POSIXdate_plotting", "cumu.DoY", "process_year", "GPP")




  rm(POSIX,  cumu.DoY, proc_yr, GPP)

  # need: process_year, GPPsat_pred, cumu.DoY
  maxGPPsat_byYear <- df %>%
    group_by(process_year) %>%
    filter(GPP == max(GPP, na.rm = TRUE)) %>%
    rename("max_GPP" = GPP) %>%
    mutate("start_cumu.DoY" = cumu.DoY - 10,
           "end_cumu.DoY" = cumu.DoY + 10)


  # need: process_year, DoY, cumu.DoY, fracyr_null, k_POSIXdate_plotting

  # This little chunk of code will calculate the maxGPPsat around +/- 10 days from the actual maxGPPsat
  maxGPPsat_byYear <- maxGPPsat_byYear %>%
    group_by(process_year) %>%
    rename("cumu.DoY_at_maxGPP" = cumu.DoY, "POSIXdate_at_maxGPP" = k_POSIXdate_plotting) %>%
    right_join(df, by = "process_year") %>%
    group_by(process_year) %>%
    subset(cumu.DoY >= start_cumu.DoY &
             cumu.DoY <= end_cumu.DoY ) %>%
    summarize("k_POSIXdate_plotting" = unique(POSIXdate_at_maxGPP),
              "cumu.DoY" = unique(cumu.DoY_at_maxGPP),
              "max_GPP_21dAvg" = mean(GPP, na.rm = TRUE),
              "max_GPP" = unique(max_GPP), .groups = "drop_last")


  SOS10_bySite <- maxGPPsat_byYear %>%
    rename("POSIXdate_at_maxGPP" = k_POSIXdate_plotting, "cumu.DoY_at_maxGPP" = cumu.DoY) %>%
    right_join(df, by = c("process_year")) %>%
    arrange(process_year, cumu.DoY) %>%
    group_by(process_year) %>%
    subset(GPP >= (0.1*max_GPP_21dAvg) &
             cumu.DoY <= cumu.DoY_at_maxGPP) %>%
    arrange(process_year, desc(cumu.DoY)) %>%
    mutate(gap = cumu.DoY - lead(cumu.DoY, default = 0)) %>%
    filter(gap != 1) %>%
    arrange(process_year, desc(cumu.DoY)) %>%
    filter(row_number()==1) %>%
    select(!c(gap, POSIXdate_at_maxGPP , cumu.DoY_at_maxGPP, max_GPP_21dAvg, max_GPP)) %>%
    mutate("SOS10" = yday(k_POSIXdate_plotting))

  SOS25_bySite <- maxGPPsat_byYear %>%
    rename("POSIXdate_at_maxGPP" = k_POSIXdate_plotting, "cumu.DoY_at_maxGPP" = cumu.DoY) %>%
    right_join(df, by = c("process_year")) %>%
    arrange(process_year, cumu.DoY) %>%
    group_by(process_year) %>%
    subset(GPP >= (0.25*max_GPP_21dAvg) &
             cumu.DoY <= cumu.DoY_at_maxGPP) %>%
    arrange(process_year, desc(cumu.DoY)) %>%
    mutate(gap = cumu.DoY - lead(cumu.DoY, default = 0)) %>%
    filter(gap != 1) %>%
    arrange(process_year, desc(cumu.DoY)) %>%
    filter(row_number()==1) %>%
    select(!c(gap, POSIXdate_at_maxGPP , cumu.DoY_at_maxGPP, max_GPP_21dAvg, max_GPP)) %>%
    mutate("SOS25" = yday(k_POSIXdate_plotting))

  SOS50_bySite <- maxGPPsat_byYear %>%
    rename("POSIXdate_at_maxGPP" = k_POSIXdate_plotting, "cumu.DoY_at_maxGPP" = cumu.DoY) %>%
    right_join(df, by = c("process_year")) %>%
    arrange(process_year, cumu.DoY) %>%
    group_by(process_year) %>%
    subset(GPP >= (0.5*max_GPP_21dAvg) &
             cumu.DoY <= cumu.DoY_at_maxGPP) %>%
    arrange(process_year, desc(cumu.DoY)) %>%
    mutate(gap = cumu.DoY - lead(cumu.DoY, default = 0)) %>%
    filter(gap != 1) %>%
    arrange(process_year, desc(cumu.DoY)) %>%
    filter(row_number()==1) %>%
    select(!c(gap, POSIXdate_at_maxGPP , cumu.DoY_at_maxGPP, max_GPP_21dAvg, max_GPP)) %>%
    mutate("SOS50" = yday(k_POSIXdate_plotting))




  EOS10_bySite <- maxGPPsat_byYear %>%
    rename("POSIXdate_at_maxGPP" = k_POSIXdate_plotting, "cumu.DoY_at_maxGPP" = cumu.DoY) %>%
    right_join(df, by = c("process_year")) %>%
    arrange(process_year, cumu.DoY) %>%
    group_by(process_year) %>%
    subset(GPP >= (0.1*max_GPP_21dAvg) &
             cumu.DoY >= cumu.DoY_at_maxGPP) %>%
    arrange(process_year, cumu.DoY) %>%
    mutate(gap = lead(cumu.DoY, default = 0) - cumu.DoY) %>%
    filter(gap != 1) %>%
    arrange(process_year, cumu.DoY) %>%
    filter(row_number()==1) %>%
    select(!c(gap, POSIXdate_at_maxGPP , cumu.DoY_at_maxGPP, max_GPP_21dAvg, max_GPP)) %>%
    mutate("EOS10" = yday(k_POSIXdate_plotting))

  EOS25_bySite <- maxGPPsat_byYear %>%
    rename("POSIXdate_at_maxGPP" = k_POSIXdate_plotting, "cumu.DoY_at_maxGPP" = cumu.DoY) %>%
    right_join(df, by = c("process_year")) %>%
    arrange(process_year, cumu.DoY) %>%
    group_by(process_year) %>%
    subset(GPP >= (0.25*max_GPP_21dAvg) &
             cumu.DoY >= cumu.DoY_at_maxGPP) %>%
    arrange(process_year, cumu.DoY) %>%
    mutate(gap = lead(cumu.DoY, default = 0) - cumu.DoY) %>%
    filter(gap != 1) %>%
    arrange(process_year, cumu.DoY) %>%
    filter(row_number()==1) %>%
    select(!c(gap, POSIXdate_at_maxGPP , cumu.DoY_at_maxGPP, max_GPP_21dAvg, max_GPP)) %>%
    mutate("EOS25" = yday(k_POSIXdate_plotting))

  EOS50_bySite <- maxGPPsat_byYear %>%
    rename("POSIXdate_at_maxGPP" = k_POSIXdate_plotting, "cumu.DoY_at_maxGPP" = cumu.DoY) %>%
    right_join(df, by = c("process_year")) %>%
    arrange(process_year, cumu.DoY) %>%
    group_by(process_year) %>%
    subset(GPP >= (0.5*max_GPP_21dAvg) &
             cumu.DoY >= cumu.DoY_at_maxGPP) %>%
    arrange(process_year, cumu.DoY) %>%
    mutate(gap = lead(cumu.DoY, default = 0) - cumu.DoY) %>%
    filter(gap != 1) %>%
    arrange(process_year, cumu.DoY) %>%
    filter(row_number()==1) %>%
    select(!c(gap, POSIXdate_at_maxGPP , cumu.DoY_at_maxGPP, max_GPP_21dAvg, max_GPP)) %>%
    mutate("EOS50" = yday(k_POSIXdate_plotting))



  CritDates_POSIX_wide <- maxGPPsat_byYear %>%
    select(process_year, k_POSIXdate_plotting, max_GPP, max_GPP_21dAvg) %>%
    rename("date_at_maxGPP" = k_POSIXdate_plotting, "maxGPP" = max_GPP, "maxGPP_21d" = max_GPP_21dAvg) %>%
    left_join(select(SOS10_bySite, process_year, k_POSIXdate_plotting, GPP), by = "process_year") %>%
    rename("date_at_SOS10" = k_POSIXdate_plotting, "GPP_at_SOS10" = GPP) %>%
    left_join(select(SOS25_bySite, process_year, k_POSIXdate_plotting, GPP), by = "process_year") %>%
    rename("date_at_SOS25" = k_POSIXdate_plotting, "GPP_at_SOS25" = GPP) %>%
    left_join(select(SOS50_bySite, process_year, k_POSIXdate_plotting, GPP), by = "process_year") %>%
    rename("date_at_SOS50" = k_POSIXdate_plotting, "GPP_at_SOS50" = GPP) %>%
    left_join(select(EOS50_bySite, process_year, k_POSIXdate_plotting, GPP), by = "process_year") %>%
    rename("date_at_EOS50" = k_POSIXdate_plotting, "GPP_at_EOS50" = GPP) %>%
    left_join(select(EOS25_bySite, process_year, k_POSIXdate_plotting, GPP), by = "process_year") %>%
    rename("date_at_EOS25" = k_POSIXdate_plotting, "GPP_at_EOS25" = GPP) %>%
    left_join(select(EOS10_bySite, process_year, k_POSIXdate_plotting, GPP), by = "process_year") %>%
    rename("date_at_EOS10" = k_POSIXdate_plotting, "GPP_at_EOS10" = GPP) %>%
    relocate(process_year, maxGPP, maxGPP_21d, date_at_SOS10, date_at_SOS25, date_at_SOS50, date_at_maxGPP, date_at_EOS50, date_at_EOS25, date_at_EOS10,
             GPP_at_SOS10, GPP_at_SOS25, GPP_at_SOS50, GPP_at_EOS50, GPP_at_EOS25, GPP_at_EOS10) %>%
    mutate(process_year = as.numeric(as.character(process_year)))


  CritDates_POSIX_long_dates <- CritDates_POSIX_wide %>%
    select(!c(GPP_at_SOS10, GPP_at_SOS25, GPP_at_SOS50, GPP_at_EOS50, GPP_at_EOS25, GPP_at_EOS10)) %>%
    rename("SOS10" = date_at_SOS10, "SOS25" = date_at_SOS25, "SOS50" = date_at_SOS50, "Peak_GPPsat" = date_at_maxGPP, "EOS50" = date_at_EOS50, "EOS25" = date_at_EOS25, "EOS10" = date_at_EOS10) %>%
    gather(key = "CritThreshold", value = "k_POSIXdate_plotting", SOS10:EOS10) %>%
    arrange(process_year, match(CritThreshold, c("SOS10","SOS25","SOS50","Peak_GPPsat","EOS50","EOS25","EOS10")))

  CritDates_POSIX_long_GPP <- CritDates_POSIX_wide %>%
    select(!c( maxGPP, date_at_SOS10, date_at_SOS25, date_at_SOS50, date_at_maxGPP, date_at_EOS50, date_at_EOS25, date_at_EOS10)) %>%
    rename("SOS10" = GPP_at_SOS10, "SOS25" = GPP_at_SOS25, "SOS50" = GPP_at_SOS50, "Peak_GPPsat" = maxGPP_21d, "EOS50" = GPP_at_EOS50, "EOS25" = GPP_at_EOS25, "EOS10" = GPP_at_EOS10) %>%
    gather(key = "CritThreshold", value = "GPP_pred",  Peak_GPPsat:EOS10) %>%
    arrange(process_year, match(CritThreshold, c("SOS10","SOS25","SOS50","Peak_GPPsat","EOS50","EOS25","EOS10")))


  CritDates_POSIX_long <- CritDates_POSIX_long_dates %>%
    left_join(CritDates_POSIX_long_GPP)

  CritDates_POSIX_long$fracyr_null <- fPOSIX_to_fracyr_null(CritDates_POSIX_long$k_POSIXdate_plotting, as.numeric(as.character(CritDates_POSIX_long$process_year)))


  CritDates <- list("CritDates_long" = CritDates_POSIX_long, "CritDates_wide" = CritDates_POSIX_wide)

  return(CritDates)
}