R/fPaginateGPPsatFacets.R

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

#' This function will plot annual facets of GPPsat time series with critical dates overlaid on each plot

#' @export
#' @title Plot faceted time series of GPPsat (with SOS/EOS)
#' @param GPP numeric, requires a vector of daily GPP for a given process_year
#' @param GPP.pred numeric, requires a vector of PREDICTED daily GPP for a given process_year
#' @param GPP.POSIX date, requires a vector of POSIX dates (YYYY-MM-DD) for each process_year
#' @param GPP.ProcYr numeric, requires a vector of process_years (used to facet our data by individual years)
#' @param crit.ID character, requires a vector of critical threshold labels (e.g. SOS10, SOS25) (vector must be of equal length as the number of individual critical dates)
#' @param crit.POSIX date, requires a vector of POSIX dates (YYYY-MM-DD) for each critical threshold
#' @param crit.ProcYr numeric, requires a vector of process_years (used to identify which subset of data the critical thresholds were estimated from)
#' @param crit.GPP numeric, requires a vector of GPPsat estimates for each critical threshold date
#' @param site.info list, metadata for current site
#' @param ncol numeric, number of columns in your faceted plot
#' @param nrow numeric, number of rows in your faceted plot
#' @param span numeric, set the span for time series smoothing (defaults to 0.075)


#' @importFrom ggforce facet_wrap_paginate



fPaginateGPPsatFacets <- function(GPP, GPP.pred, GPP.POSIX, GPP.ProcYr, crit.ID, crit.POSIX, crit.ProcYr, crit.GPP, site.info, ncol = 4, nrow = 4, span) {

  GPP.df <- data.frame(GPP.ProcYr, GPP.POSIX, GPP, GPP.pred)
  names(GPP.df) <- c("process_year", "k_POSIXdate_plotting", "GPP", "GPP_pred")


  crit.df <- data.frame(crit.ProcYr, crit.POSIX, crit.ID, crit.GPP)
  names(crit.df) <- c("process_year", "k_POSIXdate_plotting", "CriticalThreshold", "Crit_pred")

  CritDates <- crit.df %>%
    left_join(GPP.df, by = c("process_year", "k_POSIXdate_plotting")) %>%
    mutate(fracyr = decimal_date(k_POSIXdate_plotting),
           fracyr_null = fPOSIX_to_fracyr_null(k_POSIXdate_plotting, process_year),
           DoY = yday(k_POSIXdate_plotting))


  PlottingData <- GPP.df %>%
    mutate(fracyr = decimal_date(k_POSIXdate_plotting),
           fracyr_null = fPOSIX_to_fracyr_null(k_POSIXdate_plotting, process_year),
           DoY = yday(k_POSIXdate_plotting))



  # rm(GPP, GPP.pred, GPP.POSIX, GPP.ProcYr, crit.ID, crit.POSIX, crit.ProcYr)




  # Specify labels and ranges for your plots
  ylab <- expression(atop("GPPsat (NEE method, fixed 'b')",paste("(µmol ", m^-2," ", s^-1, ")")))
  yrange <- range(PlottingData$GPP, na.rm = TRUE)
  ymax <- yrange[2] + (yrange[2] - yrange[1])/10

  # Set a limits to the yaxis range that's reasonable (if your data falls below this threshold, it won't be plotted)
  if (yrange[1] < -10) {
    yrange[1] <- -10
  }

  if (ymax > 40) {
    ymax <- 40
  }

  ylim <- c(yrange[1], ymax)

  xmin <- min(PlottingData$fracyr_null)
  xmax <- max(PlottingData$fracyr_null)


  colors <- c("ECW" = "#e7d4e8", "SOS10" = "#1a9850", "SOS25" = "#91cf60", "SOS50" = "#d9ef8b", "Peak_GPPsat" = "white",
              "EOS10" = "#d73027", "EOS25" = "#fc8d59", "EOS50" = "#fee08b", "SNW" = "#762a83")

  colors <- colors[order(factor(names(colors), levels = c("ECW", "SOS10", "SOS25", "SOS50", "Peak_GPPsat", "EOS50", "EOS25", "EOS10", "SNW")))]


  # Specify breaks for the x-axis
  DoYBreaks <- fracyr_to_DoYBreaks(interval = 100, xmin = xmin, xmax = xmax, padding = FALSE )



  plot_base <- ggplot(data = PlottingData) +
    geom_hline(yintercept = 0, color = "gray") +
    geom_point(aes(x = fracyr_null, y = GPP), alpha = 0.1) +
    geom_line(aes(x = fracyr_null, y = GPP_pred), size = 0.5, alpha = 0.5) +
    scale_y_continuous(name= ylab, limits = ylim) +
    scale_x_continuous(name = "DoY", breaks = DoYBreaks$fracyrBreaks,
                       labels = DoYBreaks$DoYBreaks,
                       limits = c(DoYBreaks$lwr[1], DoYBreaks$upr[1])) +
    ggtitle(paste0(site.info$site, ", ", site.info$db, " (", site.info$info, ")"),
            subtitle = paste0("Country: ", site.info$country, "          lat: ", round(site.info$lat, 2),"; long: ", round(site.info$long,2),
                              "          MAT: ", round(site.info$MAT, 1),"°C; MAP: ", round(site.info$MAP,0), "mm",
                              "          Forest type: ", site.info$IGBP, "\n\nspan = ", span))  +
    theme_time() +
    annotate("segment", x=-Inf, xend=Inf, y=-Inf, yend=-Inf, size=1.5, color = "black") +
    annotate("segment", x=-Inf, xend=-Inf, y=-Inf, yend=Inf, size=1.5, color = "black") +
    theme(plot.title = element_text(color="black", hjust = 0.5, size=14, face="bold"),
          plot.subtitle = element_text(color="black", hjust = 0.5, size=12),
          strip.background = element_blank(),
          strip.text = element_text(size = 12, face = "bold"),
          axis.line = element_blank(),
          legend.position = "bottom",
          legend.title = element_blank(),
          panel.grid.major = element_line(linetype = "dotted", color = "gray")) +
    guides(fill = guide_legend(nrow = 1))


  plot_base <- plot_base +
    geom_point(data = CritDates, aes(x = fracyr_null, y = GPP_pred, fill = CriticalThreshold), shape = 21, color = "black", size = 3) +
    geom_text(data = CritDates, aes(x = fracyr_null, y = GPP_pred, label = DoY),
                nudge_x = 0.15, vjust = 0, size = 3) +
    scale_fill_manual(name = "Critical Date", values = colors,
                      breaks = c("SOS10", "SOS25", "SOS50","Peak_GPPsat", "EOS50", "EOS25", "EOS10"),
                      labels = c(expression(SOS[10]), expression(SOS[25]), expression(SOS[50]),
                                 expression(maxGPP["sat"]), expression(EOS[50]), expression(EOS[25]), expression(EOS[10]) ))




  # if you have 3 x 3 facets, then you have 9 facets per page.
  # If you have 20 years of data, that's 3 pages: two have 9 facets, the third has 2 facets


  no_dat_years <- length(unique(PlottingData$process_year))
  no_facets <- nrow * ncol
  no_pages <- ceiling(no_dat_years/no_facets)

  plot_facet_list <- list()


  # for (i in 1:no_pages) {
  #
  #   plot_facet_i <- eval(substitute(
  #     plot_base +
  #     facet_wrap_paginate(. ~ process_year, ncol = ncol, nrow = ncol, page = i) ))
  #
  #
  #   print(i)
  #   print(plot_facet_i)
  #   plot_facet_list[[i]] <- plot_facet_i
  #
  # }

  for (i in 1:no_pages) {

    plot_facet_i <- plot_base +
      facet_wrap_paginate(. ~ process_year, ncol = ncol, nrow = nrow, page = i)

    plot_facet_list[[i]] <- ggplotGrob(plot_facet_i)

  }



  class(plot_facet_list) <- c("arrangelist", class(plot_facet_list))

  return(plot_facet_list)


}