FluxSynthU: Suite of functions to analyze eddy covariance data

Documented in fPlot_FullTimeSeries_16panel_InitialQC

#' This function takes a dataframe of timestamped environmental/flux variables and plots them as individual facets in a 16-panel plot.
#' Data are plotted as a continuous time series from the start of the record through the end (with no time gaps).


#' @export
#' @title 16-panel composite plot for initial data QC (continuous time series)
#' @param site.info a named vector containing the following variables: site, db, info, country, lat, long, MAT, MAP, IGBP
#' @param dat a dataframe containing timestamped variable columns
#' @param POSIXdate_col a character vector containing POSIX dates in the format 'YYYY-MM-DD' (can be either class 'date' or 'character')
#'
#' @details A bit more information about the 'site.info' dataframe columns:
#' \describe{
#' \item{site}{character object denoting the abbreviated site name (e.g. 'US-NR1')}
#' \item{db}{character object denoting the abbreviated database name (e.g. 'FLX')}
#' \item{info}{character string denoting the full site name (e.g. 'Niwot Ridge, CO')}
#' \item{country}{character object denoting the country of origin for site (e.g. 'United States')}
#' \item{lat}{numeric object denoting the site latitude}
#' \item{long}{numeric object denoting the site longitude}
#' \item{MAT}{numeric object denoting the site mean annual temperature (deg C)}
#' \item{MAP}{numeric object denoting the site mean annual precipitation (mm)}
#' \item{IGBP}{character object denoting the International Geosphere-Biosphere Programme land classification (e.g. 'ENF')}
#' }
#' @importFrom ggplot2 ggplot aes geom_line geom_vline scale_x_date scale_y_continuous ggtitle theme element_text element_blank element_line



fPlot_FullTimeSeries_16panel_InitialQC <- function(site.info, dat, POSIXdate_col){

  # Function to determine whether POSIXdate_col is of class 'Date'
  is.Date <- function(x) inherits(x, "Date")

  # Make sure your 'x' variable is of 'Date' class and matches the date format 'YYYY-MM-DD'
  if (!is.Date(dat[,POSIXdate_col])) {
    message("Converting 'POSIXdate_col' to class 'Date'")
    dat[,POSIXdate_col] <- as.Date(dat[,POSIXdate_col])

  }

  # Create your limits and breaks right here
  # Plotting limits (here we're adding 1 month to either end of the record)
  plotting_start_date <- seq(dat[[POSIXdate_col]][1], length = 2, by = "-1 months")[2]
  plotting_end_date <- seq(tail(dat[[POSIXdate_col]],1), length = 2, by = "+1 months")[2]
  date.lims <-  c(plotting_start_date, plotting_end_date)

  years_of_record <- seq(from = lubridate::year(date.lims[1]), to = lubridate::year(date.lims[2]), by = 1)

  date.range <- as.Date(c(paste(years_of_record[1], 1, 1, sep = "-"),
                          paste(year(plotting_end_date), 1, 1, sep = "-")))


  # Need 3 things:
  # breaks
  # break labels
  # label format


  if (length(years_of_record) < 3) {
    month.breaks <-  seq.Date(date.range[1], date.range[2], by = "month")
    year.breaks <- seq.Date(date.range[1], date.range[2], by = "year")

    date.breaks <- seq.Date(date.range[1], date.range[2], by = "4 months")

    date.labs <- format(date.breaks, "%Y-%m")

  } else if (length(years_of_record) < 6) {
    year.breaks <- seq.Date(date.range[1], date.range[2], by = "year")
    date.breaks <- year.breaks

    date.labs <- format(date.breaks, "%Y")

  } else if (length(years_of_record) < 9) {
    year.breaks <- seq.Date(date.range[1], date.range[2], by = "year")
    date.breaks <- year.breaks

    date.labs <- format(date.breaks, "'%y")

  } else {
    year.breaks <- seq.Date(date.range[1], date.range[2], by = "year")
    date.breaks <- year.breaks[c(TRUE,FALSE)]

    date.labs <- format(date.breaks, "'%y")
  }


  # plot time series of: k_NEE
  f.NEE_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_NEE") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("NEE",paste("(µmol ",CO[2]," ",m^-2," ",s^-1,")")))) +
    ggtitle("Net Ecosystem Exchange")

  # plot time series of: k_LE
  f.LE_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_LE") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name = expression(atop("LE", paste("(W ",m^-2,")")))) +
    ggtitle("Latent Heat Flux")

  # plot time series of: k_H
  f.H_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_H") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("H", paste("(W ",m^-2,")")))) +
    ggtitle("Sensible Heat Flux")

  # plot time series of: k_SW_in
  f.SW_in_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_SW_in") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("SW_in", paste("(W ",m^-2,")")))) +
    ggtitle("SW radiation, incoming")

  # plot time series of: k_SW_out
  f.SW_out_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_SW_out") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("SW_out", paste("(W ",m^-2,")")))) +
    ggtitle("SW radiation, outgoing")

  # plot time series of: k_LW_in
  f.LW_in_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_LW_in") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("LW_in", paste("(W ",m^-2,")")))) +
    ggtitle("LW radiation, incoming")

  # plot time series of: k_LW_out
  f.LW_out_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_LW_out") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("LW_out", paste("(W ",m^-2,")")))) +
    ggtitle("LW radiation, outgoing")

  # plot time series of: k_albedo
  f.albedo_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_albedo") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("Albedo")), limits = c(0,1)) +
    ggtitle("Albedo")

  # plot time series of: k_Rnet
  f.Rnet_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_Rnet") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("Rnet", paste("(W ",m^-2,")")))) +
    ggtitle("Net radiation")

  # plot time series of: k_Tair
  f.Tair_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_Tair") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("Tair", paste("(\u00B0C)")))) +
    ggtitle("Air temperature")

  # plot time series of: k_Tsoil
  f.Tsoil_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_Tsoil") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("Tsoil", paste("(\u00B0C)")))) +
    ggtitle("Soil temperature")

  # plot time series of: k_SWC
  f.SWC_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_SWC") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("SWC", paste("(%)")))) +
    ggtitle("Soil water content")

  # plot time series of: k_RH
  f.RH_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_RH") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("RH", paste("(%)")))) +
    ggtitle("Relative humidity")

  # plot time series of: k_VPD
  f.VPD_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_VPD") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("VPD", paste("(hPa)")))) +
    ggtitle("Vapor pressure deficit")

  # plot time series of: k_ustar
  f.ustar_full <-fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_ustar") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("U*", paste("(",m^-2,")")))) +
    ggtitle("Friction velocity")

  #plot time series of: k_PPFD
  f.PPFD_full <- fPlot_Line_Single(dat = dat, xcol = POSIXdate_col, ycol = "k_PPFD") +
    geom_vline(xintercept=year.breaks,
               linetype="longdash", alpha=0.4, color="black", size = 0.2) +
    scale_x_date(breaks = date.breaks,
                 limits = date.lims,
                 date_labels = date.labs) +
    scale_y_continuous(name=expression(atop("PPFD",paste("(µmol photon ",m^-2," ",s^-1,")")))) +
    ggtitle("Photosynthetic photon flux density")

  # Create faceted plot of all subplots
  f.composite_full <- f.NEE_full + f.LE_full + f.H_full + f.Rnet_full + f.SW_in_full + f.SW_out_full + f.LW_in_full + f.LW_out_full + f.albedo_full +
    f.Tair_full + f.Tsoil_full + f.SWC_full + f.RH_full + f.VPD_full + f.ustar_full + f.PPFD_full + plot_layout(ncol=4)


  fig_full <- patchwork::wrap_elements(f.composite_full)

  # Add a title to your plot that includes all of the site metadata
  fig_full <- fAddSiteMeta2Plot(fig_full, site = site.info$site,
                                db = site.info$db,
                                longname = site.info$info,
                                country = site.info$country,
                                lat = site.info$lat,
                                long = site.info$long,
                                MAT = site.info$MAT,
                                MAP = site.info$MAP,
                                IGBP = site.info$IGBP)

  newlist <- list(fig_full, f.NEE_full, f.LE_full, f.H_full, f.Rnet_full, f.SW_in_full, f.SW_out_full, f.LW_in_full, f.LW_out_full, f.albedo_full,
                  f.Tair_full, f.Tsoil_full, f.SWC_full, f.RH_full, f.VPD_full, f.ustar_full, f.PPFD_full)
  names(newlist) <- c("fig_full", "f.NEE_full", "f.LE_full", "f.H_full", "f.Rnet_full", "f.SW_in_full", "f.SW_out_full", "f.LW_in_full", "f.LW_out_full", "f.albedo_full",
                      "f.Tair_full", "f.Tsoil_full", "f.SWC_full", "f.RH_full", "f.VPD_full", "f.ustar_full", "f.PPFD_full")
  return(newlist)
}