R/visualization.R

In grahamstewart12/tidyflux:

## VISUALIZATION

plot_rolling_corrs <- function(data, x, vars = c("tair", "vpd", "ppfd"),
                               n = 336, remove_outliers = TRUE, z = 7,
                               facets = FALSE, show_vars = FALSE) {
  if (isTRUE(show_vars) & isFALSE(facets)) {
    stop("'show_vars' cannot be used without 'facets'.")
  }
  diff <- as.numeric(diff(data$timestamp)[1])
  # Assuming time difference is 1 hour if it is not 30 mins
  days <- if (diff == 30) n / 48 else n / 24
  if (days < 7) warning(
    "Rolling correlations probably not meaningful on sub-weekly time periods."
  )
  x <- rlang::enquo(x)
  name <- data %>% dplyr::select(!!x) %>% names()
  data <- data %>%
    dplyr::select(timestamp, !!x, vars) %>%
    dplyr::mutate(
      ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA),
      x_hack = !!x # complete hack to issue w/ unquoting in roll_corr()
    ) %>%
    dplyr::mutate_at(
      vars, ~roll_corr(x_hack, ., n)
    ) %>%
    #dplyr::mutate_at(vars, .funs = list(cat = ~ntile(., 2)))
    tidyr::gather(key = "key", value = "value", -timestamp, -!!x, -x_hack) %>%
    dplyr::rowwise() %>%
    dplyr::mutate(key = get_full_names()[[key]]) %>%
    dplyr::ungroup() %>%
    dplyr::group_by(key) %>%
    dplyr::mutate(
      # Get rid of first and last corr values - they are usually bad
      value = replace(value, seq_along(value) %in% c(
        which(!is.na(value))[1], rev(which(!is.na(value)))[1]
      ), NA),
      sign = dplyr::if_else(value >= 0, "pos", "neg")
    ) %>%
    dplyr::ungroup()
  # Create caption
  cap <- paste0("Correlations calculated on a rolling basis with n = ", n)
  # Create base plot
  out <- data %>%
    ggplot2::ggplot(ggplot2::aes(timestamp, value)) +
    ggplot2::geom_hline(aes(yintercept = 0), size = 0.25, color = "grey60") +
    ggplot2::theme_bw() +
    ggplot2::labs(
      y = paste0("Pearson correlation (with ", get_full_names()[[name]], ")"),
      x = NULL
    )
  # Add geoms based on whether 'facets' is indicated
  if (facets) {
    out <- out +
      geom_line(aes(color = value), na.rm = TRUE, size = 0.75) +
      scale_color_gradient2(
        low = "steelblue", mid = "grey90", high = "orangered"
      ) +
      #ggplot2::geom_bar(
      #  ggplot2::aes(fill = sign), stat = "identity", na.rm = TRUE
      #) +
      #ggplot2::scale_fill_manual(
      #  values = c("steelblue", "orangered"),
      #  labels = c("night", "day")
      #) +
      ggplot2::theme(legend.position = "none") +
      ggplot2::facet_wrap(~key, ncol = 1)
  } else {
    out <- out + ggplot2::geom_line(aes(color = key), na.rm = TRUE)
  }
  out
}

plot_correlations <- function(data, omit, remove_outliers = TRUE, z = 7,
                              show_coefs = FALSE) {
  if (!missing(omit)) omit <- rlang::enquos(omit)
  vars <- c(
    "xpeak", "swc", "wd", "pair", "fh2o", "vpd", "ppfd", "h", "ws", "fco2",
    "fch4", "tair", "tsoil"
  )
  #data <- data %>% dplyr::select(-timestamp, -daytime)
  data <- data %>% dplyr::select(vars)
  cormat <- round(cor(data, use = "pairwise.complete.obs"), 2)
  # Get lower triangle of the correlation matrix
  get_lower_tri <- function(x) {
    x[upper.tri(x, diag = TRUE)] <- NA
    return(x)
  }
  # Get upper triangle of the correlation matrix
  get_upper_tri <- function(x){
    x[lower.tri(x, diag = TRUE)] <- NA
    return(x)
  }
  # Helper function to reorder the correlation matrix
  reorder_cormat <- function(x) {
    # Use correlation between variables as distance
    dd <- as.dist((1 - x) / 2)
    hc <- hclust(dd)
    x <- x[hc$order, hc$order]
  }
  # Reorder the correlation matrix
  cormat <- reorder_cormat(cormat)
  upper_tri <- get_upper_tri(cormat)
  # Melt the correlation matrix
  melted_cormat <- reshape2::melt(upper_tri, na.rm = TRUE)
  # Create a ggheatmap
  out <- melted_cormat %>%
    ggplot2::ggplot(ggplot2::aes(Var2, Var1, fill = value)) +
    ggplot2::geom_tile(color = "grey80", size = 0.3) +
    ggplot2::scale_fill_gradient2(
      low = "steelblue", high = "orangered", midpoint = 0, limit = c(-1, 1),
      space = "Lab", name = "Pearson\nCorrelation"
    ) +
    ggplot2::labs(x = NULL, y = NULL) +
    ggplot2::theme_minimal() +
    ggplot2::theme(
      legend.position = c(0, 1),
      legend.justification = c(0, 1),
      axis.text.x = ggplot2::element_text(
        angle = 45, vjust = 1, size = 12, hjust = 1
      )
    ) +
    ggplot2::coord_fixed()
  # Add correlation coefficients
  if (show_coefs) {
    out <- out +
      ggplot2::geom_text(
        ggplot2::aes(Var2, Var1, label = value), color = "black", size = 4
      ) +
      ggplot2::theme(
        axis.title.x = ggplot2::element_blank(),
        axis.title.y = ggplot2::element_blank(),
        panel.grid.major = ggplot2::element_blank(),
        panel.border = ggplot2::element_blank(),
        panel.background = ggplot2::element_blank(),
        axis.ticks = ggplot2::element_blank(),
        legend.justification = c(1, 0),
        legend.position = c(0.6, 0.7),
        legend.direction = "horizontal") +
      ggplot2::guides(fill = ggplot2::guide_colorbar(
        barwidth = 7, barheight = 1, title.position = "top", title.hjust = 0.5)
      )
  }
  out
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param remove_outliers
#' @param z
#' @param groups
#'
#' @return
#' @export
#'
#' @examples
plot_distribution <- function(data, x, remove_outliers = TRUE, z = 7,
                              groups = c("none", "daynight")) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()

  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n())
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n() - n_all) %>%
      unlist()
  }
  if (groups == "daynight") {
    data <- data %>% dplyr::mutate(daytime = factor(daytime))
  }

  # Build plot caption according to how data were processed
  caption <- NULL
  if (remove_outliers) {
    caption <- paste0(n_removed, " outliers removed using MAD with z = ", z)
  }

  # Plot
  data %>%
    tidyr::drop_na() %>%
    ggplot2::ggplot(ggplot2::aes(!!x)) +
    {if (groups == "daynight") {
      ggplot2::geom_density(
        ggplot2::aes(fill = daytime, color = daytime), alpha = 0.3
      )
    }} +
    {if (groups == "none") {
      ggplot2::geom_density(alpha = 0.3, fill = "grey80")
    }} +
    {if (groups == "daynight") {
      ggplot2::scale_color_manual(
        values = c("steelblue", "orangered"),
        labels = c("night", "day")
      )
    }} +
    {if (groups == "daynight") {
      ggplot2::scale_fill_manual(
        values = c("steelblue", "orangered"),
        labels = c("night", "day")
      )
    }} +
    ggplot2::labs(
      y = NULL, x = paste0(name, " (", units, ")"), caption = caption
    ) +
    ggplot2::theme_bw() +
    ggplot2::theme(
      axis.text.y = ggplot2::element_blank(),
      axis.ticks.y = ggplot2::element_blank(),
      legend.position = c(1, 1),
      legend.justification = c(1, 1),
      legend.key = ggplot2::element_rect(
        fill = "transparent", colour = "transparent"
      ),
      legend.background = ggplot2::element_rect(
        fill = "transparent", colour = "transparent"
      ),
      legend.title = ggplot2::element_blank(),
      legend.spacing.x = unit(0.1, "cm"),
      #legend.box = "horizontal",
      legend.key.size = unit(1, "line")
    )
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param timestamp
#' @param tair
#' @param tsoil
#' @param p
#' @param remove_outliers
#' @param z
#' @param smooth
#' @param n
#' @param title
#'
#' @return
#' @export
#'
#' @examples
plot_daily_means <- function(data, x, timestamp = timestamp, tair = tair,
                             tsoil = tsoil, p = p, remove_outliers = TRUE,
                             z = 7, smooth = TRUE, n = 3, title = FALSE) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  timestamp <- rlang::enquo(timestamp)
  tair <- rlang::enquo(tair)
  tsoil <- rlang::enquo(tsoil)
  p <- rlang::enquo(p)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()
  t_units <- data %>% dplyr::pull(!!tair) %>% attr("units")
  p_units <- data %>% dplyr::pull(!!p) %>% attr("units")

  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n())
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n() - n_all) %>%
      unlist()
  }

  # Extract time components from timestamp
  data <- data %>%
    dplyr::mutate(date = lubridate::date(!!timestamp)) %>%
    dplyr::mutate(month = lubridate::month(!!timestamp, label = TRUE))

  # Create daily-grouped dataset, get mean
  daily <- data %>%
    dplyr::select(date, daytime, !!x, !!tair, !!tsoil, !!p) %>%
    # Extract day/night subsets of x
    dplyr::mutate(
      x_day = !!x,
      x_day = replace(x_day, is.na(daytime) | daytime == 0, NA),
      x_night = !!x,
      x_night = replace(x_night, is.na(daytime) | daytime == 1, NA)
    ) %>%
    dplyr::mutate(
      !!rlang::as_name(p) := !!p * 500
    ) %>%
    dplyr::group_by(date) %>%
    dplyr::summarize(
      # Need at least 6 values for mean (else NA)
      n_day = sum(!is.na(x_day)),
      n_night = sum(!is.na(x_night)),
      nrow = sum(!is.na(!!x)),
      x_day = dplyr::if_else(n_day > 5, mean(x_day, na.rm = TRUE), NA_real_),
      x_night = dplyr::if_else(
        n_night > 5, mean(x_night, na.rm = TRUE), NA_real_
      ),
      !!rlang::as_name(x) := dplyr::if_else(
        nrow > 5, mean(!!x, na.rm = TRUE), NA_real_
      ),
      !!rlang::as_name(tair) := mean(!!tair, na.rm = TRUE),
      !!rlang::as_name(tsoil) := mean(!!tsoil, na.rm = TRUE),
      !!rlang::as_name(p) := sum(!!p, na.rm = TRUE)
    ) %>%
    # Smooth using moving average (if requested)
    {if (smooth) {
      dplyr::mutate(
        ., x_day = zoo::rollapply(x_day, n, mean, na.rm = TRUE, fill = NA),
        x_night = zoo::rollapply(x_night, n, mean, na.rm = TRUE, fill = NA),
        !!rlang::as_name(x) := zoo::rollapply(
          !!x, n, mean, na.rm = TRUE, fill = NA
        ),
        !!rlang::as_name(tair) := zoo::rollapply(
          !!tair, n, mean, na.rm = TRUE, fill = NA
        ),
        !!rlang::as_name(tsoil) := zoo::rollapply(
          !!tsoil, n, mean, na.rm = TRUE, fill = NA
        )
      )
    } else .} %>%
    # Shift dates to improve plotting
    dplyr::mutate(date = lubridate::parse_date_time(date, "Ymd") + 43200)

  # Calculate y-axis scaling factor for p
  max <- daily %>%
    dplyr::ungroup() %>%
    dplyr::summarize(max(!!p, na.rm = TRUE)) %>%
    unlist()

  # Build plot caption according to how data were processed
  caption <- NULL
  if (remove_outliers) {
    caption <- paste0(n_removed, " outliers removed using MAD with z = ", z)
  }
  if (smooth) {
    caption <- paste0(
      caption, "; daily means smoothed using moving average with n = ", n
    )
  }

  # Create time series plot
  flux <- daily %>%
    ggplot2::ggplot(ggplot2::aes(date, !!x)) +
    ggplot2::geom_line(
      ggplot2::aes(y = x_day, linetype = "day"), na.rm = TRUE, alpha = 0.5
    ) +
    ggplot2::geom_line(
      ggplot2::aes(y = x_night, linetype = "night"), na.rm = TRUE, alpha = 0.5
    ) +
    ggplot2::geom_line(ggplot2::aes(y = !!x, linetype = "all"), na.rm = TRUE) +
    ggplot2::labs(
      x = NULL, y = paste0(name, " (", units, ")")
    ) +
    ggplot2::scale_linetype_manual(
      values = c("all" = 1, "day" = 2, "night" = 3),
      labels = c("all", "day", "night")
    ) +
    # Add caption
    ggplot2::labs(caption = caption) +
    ggplot2::guides(
      linetype = ggplot2::guide_legend(ncol = 3)
    ) +
    ggplot2::theme_bw() +
    ggplot2::theme(
      legend.position = c(0.99, 0.99),
      legend.justification = c(0.99, 0.99),
      legend.margin = ggplot2::margin(0, 0, 0, 0),
      legend.box.margin = ggplot2::margin(0, 0, 0, 0),
      legend.title = ggplot2::element_blank(),
      legend.spacing.x = unit(0.1, "cm"),
      legend.box = "horizontal",
      legend.key.size = unit(1, "line")
    )

  # Create tair, tsoil, and p auxilliary plot
  met <- daily %>%
    ggplot2::ggplot(ggplot2::aes(x = date)) +
    ggplot2::geom_line(
      ggplot2::aes(y = !!tair, color = "tair"), na.rm = TRUE
    ) +
    ggplot2::geom_line(
      ggplot2::aes(y = !!tsoil, color = "tsoil"), na.rm = TRUE
    ) +
    ggplot2::geom_tile(ggplot2::aes(
      y = -1 * (!!p / 2 - max), height = !!p, fill = "p"
    ), alpha = 0.5) +
    ggplot2::geom_text(
      data = daily %>% dplyr::filter(!!p / 0.5 > 50),
      ggplot2::aes(y = -1 * (!!p / 1 - max), label = as.character(!!p / 0.5)),
      hjust = 0, vjust = 0, nudge_x = 200000, nudge_y = -3, size = 3,
      color = "cornflowerblue") +
    ggplot2::labs(x = NULL, y = paste0("temperature (", t_units, ")")) +
    ggplot2::scale_color_manual(
      values = c("tair" = "steelblue", "tsoil" = "orangered"),
      labels = c("air", "soil")
    ) +
    ggplot2::scale_fill_manual(
      values = c("p" = "cornflowerblue"),
      labels = c("precip (mm)")
    ) +
    ggplot2::guides(
      col = ggplot2::guide_legend(ncol = 2),
      fill = guide_legend(override.aes = list(size = 1))
    ) +
    {if (title) {
      ggplot2::labs(subtitle = "Smoothed daily means with weather conditions")
    }} +
    ggplot2::theme_bw() +
    ggplot2::theme(
      axis.text.x = ggplot2::element_blank(),
      axis.ticks.x = ggplot2::element_blank(),
      legend.position = c(0.99, 0.99),
      legend.justification = c(0.99, 0.99),
      #legend.key = ggplot2::element_rect(
      #  fill = "transparent", colour = "transparent"
      #),
      #legend.background = ggplot2::element_rect(
      #  fill = "transparent", colour = "transparent"
      #),
      legend.margin = ggplot2::margin(0, 0, 0, 0),
      legend.box.margin = ggplot2::margin(0, 0, 0, 0),
      legend.title = ggplot2::element_blank(),
      legend.spacing.x = unit(0.1, "cm"),
      legend.box = "horizontal",
      legend.key.size = unit(1, "line")
    )

  # Return both plots stacked
  cowplot::plot_grid(
    met, flux, nrow = 2, align = "v", rel_heights = c(2/5, 3/5)
  )
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param max_covs
#' @param remove_outliers
#' @param z
#' @param filter_qc
#' @param qc_var
#' @param max_qc
#' @param show_regression
#' @param formula
#' @param sample
#' @param perc
#' @param pred_vars
#'
#' @return
#' @export
#'
#' @examples
plot_covariates <- function(data, x, max_covs = 4, remove_outliers = TRUE,
                            z = 7,  filter_qc = FALSE, qc_var, max_qc = 0,
                            show_regression = TRUE, formula = y ~ poly(x, 2),
                            sample = FALSE, perc = 0.25,
                            pred_vars = c(vpd, swc, tsoil, g, tair, ppfd, rh,
                                          ws, wd, h, le, tau, fco2, fch4, pair)
                            ) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()
  all_units <- tidyflux::units(data, "names")

  # Filter by quality control flag (if requested)
  if (filter_qc) {
    if (missing(qc_var)) {
      qc_var <- paste0("qc_", name)
    } else qc_var <- rlang::ensym(qc_var)
    # Calculate initial number of points
    n_all <- data %>%
      dplyr::filter(!is.na(!!x)) %>%
      dplyr::summarize(dplyr::n()) %>%
      unlist()
    data <- data %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, !!as.name(qc_var) > max_qc, NA)
      )
    # Calculate number of points remaining
    n_left <- data %>%
      dplyr::filter(!is.na(!!x)) %>%
      dplyr::summarize(dplyr::n()) %>%
      unlist()
  }

  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all2 <- data %>%
      dplyr::filter(!is.na(!!x)) %>%
      dplyr::summarize(dplyr::n()) %>%
      unlist()
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(!is.na(!!x)) %>%
      dplyr::summarize(n_all2 - dplyr::n()) %>%
      unlist()
  }

  # Get significant covariates
  vars <- data %>%
    dplyr::select(
      !!x, vpd, swc, tsoil, g, rn, tair, ppfd, rh, ws, wd, pair
    ) %>%
    {suppressWarnings(tidyr::gather(., "var", "value", -!!x))} %>%
    tidyr::nest(., -var) %>%
    dplyr::mutate(
      data = purrr::map(data, ~ dplyr::filter(.,
        !get_outliers(!!x, .x$value)
      )),
      test = purrr::map(
        data, ~ cor.test(~ !!x + value, data = .x, na.action = na.exclude)
      ),
      tidied = purrr::map(test, broom::tidy)
    ) %>%
    tidyr::unnest(tidied, .drop = TRUE) %>%
    dplyr::arrange(p.value) %>%
    dplyr::slice(1:max_covs) %>%
    dplyr::pull(var)

  # Plot
  data %>%
    dplyr::select(c(name, vars)) %>%
    {suppressWarnings(tidyr::gather(., "var", "value", -!!x))} %>%
    dplyr::rowwise() %>%
    dplyr::mutate(var = paste0(
      get_full_names()[[var]], " (", tidyflux::units(data[[var]]), ")"
    )) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(var = factor(var)) %>%
    dplyr::filter(!is.na(value) & !is.na(!!x)) %>%
    {if (sample) {
      dplyr::group_by(., var)
    } else .} %>%
    {if (sample) {
      dplyr::sample_frac(., perc)
    } else .} %>%
    {if (sample) {
      dplyr::ungroup(.)
    } else .} %>%
    ggplot2::ggplot(ggplot2::aes(value, !!x)) +
    ggplot2::geom_point(na.rm = TRUE, alpha = 0.3, size = 0.75) +
    ggplot2::geom_smooth(
      method = "lm", formula = formula, color = "orangered",
      fill = "steelblue", alpha = 0.5, se = FALSE
    ) +
    ggplot2::labs(x = NULL, y = paste0(name, " (", units, ")")) +
    # Add caption if outliers are removed
    {if (remove_outliers) {
      ggplot2::labs(
        caption = paste0(n_removed, " outliers removed using MAD with z = ", z)
      )
    }} +
    # Add caption if quality filter is applied
    {if (filter_qc) {
      ggplot2::labs(caption = paste0(
        "Showing ", n_left, " of ", n_all, " values: ",
        "data filtered by ", qc_var, " <= ", max_qc
      ))
    }} +
    ggplot2::facet_wrap(
      ~var, scales = "free_x", nrow = 1, labeller = ggplot2::label_wrap_gen()
    ) +
    ggplot2::theme_bw()
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param wd
#' @param bin_size
#' @param remove_outliers
#' @param z
#' @param groups
#' @param polar_coords
#' @param sample
#' @param perc
#'
#' @return
#' @export
#'
#' @examples
plot_wd_stats <- function(data, x, wd = wd, bin_size = 10,
                          remove_outliers = TRUE, z = 7,
                          groups = c("none", "daynight"),
                          polar_coords = FALSE, sample = FALSE, perc = 0.25) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  wd <- rlang::enquo(wd)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()
  groups <- match.arg(groups)
  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n())
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n() - n_all) %>%
      unlist()
  }

  # Pre-process data for plotting
  data <- data %>%
    dplyr::mutate(wd_bin = bin_wd(!!wd, bin_size = bin_size)) %>%
    # Group by daynight if indicated
    {if (groups == "daynight") {
      dplyr::filter(., !is.na(daytime)) %>%
      dplyr::mutate(
        ., daytime = recode(daytime, `0` = "Night", `1` = "Day")
      ) %>%
      # Add 3rd level to 'daytime' which includes all data
      dplyr::bind_rows(., transform(., daytime = "All")) %>%
      dplyr::group_by(., daytime)
    } else .} %>%
    {if (sample) {
      dplyr::sample_frac(., perc)
    } else .}

  # Create binned dataset
  wd_df <- data %>%
    {if (groups == "daynight") {
      dplyr::filter(., !is.na(daytime)) %>%
      # Add 3rd level to 'daytime' which includes all data
      dplyr::bind_rows(., transform(., daytime = "All")) %>%
      dplyr::group_by(., daytime, wd_bin)
    } else dplyr::group_by(., wd_bin)} %>%
    dplyr::summarize(
      mean = mean(!!x, na.rm = TRUE), sd = sd(!!x, na.rm = TRUE)
    )

  # Build plot caption according to how data were processed
  caption <- "dashed line is overall mean"
  if (remove_outliers) {
    caption = paste0(
      caption, "; ", n_removed, " outliers removed using MAD with z = ", z
    )
  }

  # Plot
  data %>%
    ggplot2::ggplot(ggplot2::aes(!!wd, !!x)) +
    # Add individual points
    ggplot2::geom_point(size = 0.75, na.rm = TRUE, alpha = 0.3) +
    # Shade area of mean +/- sd
    ggplot2::geom_ribbon(
      data = wd_df,
      ggplot2::aes(x = wd_bin, ymin = mean - sd, ymax = mean + sd),
      fill = "steelblue", alpha = 0.5, inherit.aes = FALSE
    ) +
    # Bottom shading border
    ggplot2::geom_line(
      data = wd_df,
      ggplot2::aes(wd_bin, mean - sd), color = "steelblue", na.rm = TRUE
    ) +
    # Top shading border
    ggplot2::geom_line(
      data = wd_df,
      ggplot2::aes(wd_bin, mean + sd), color = "steelblue", na.rm = TRUE
    ) +
    # Add means as points connected by lines
    ggplot2::geom_point(
      data = wd_df, ggplot2::aes(wd_bin, mean),
      color = "orangered", na.rm = TRUE
    ) +
    ggplot2::geom_line(
      data = wd_df, ggplot2::aes(wd_bin, mean),
      color = "orangered", na.rm = TRUE
    ) +
    # Add overall mean as line
    ggplot2::geom_hline(
      data = data %>%
      {if (groups == "daynight") {
        dplyr::group_by(., daytime)
      } else .} %>%
        dplyr::summarize(mean = mean(!!x, na.rm = TRUE)),
      ggplot2::aes(yintercept = mean),
      color = "orangered", linetype = 2
    ) +
    # Label x-axis with degrees
    ggplot2::scale_x_continuous(
      breaks = c(seq(0, 360, by = 100)),
      labels = c("0°", "100°", "200°", "300°")
    ) +
    ggplot2::labs(x = NULL, y = paste0(name, " (", units, ")")) +
    ggplot2::labs(caption = caption) +
    {if (polar_coords) {
      ggplot2::coord_polar()
    }} +
    {if (groups == "daynight") {
      ggplot2::facet_wrap(~daytime)
    }} +
    ggplot2::theme_bw()
}

qplot_wd_stats <- function(data, x, wd, bin_size = 10, remove_outliers = TRUE,
                           z = 7) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  wd <- rlang::enquo(wd)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()

  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
  }

  # Plot
  data %>%
    dplyr::mutate(wd_bin = bin_wd(!!wd, bin_size = bin_size)) %>%
    ggplot2::ggplot(ggplot2::aes(!!wd, !!x)) +
    # Add individual points
    ggplot2::geom_point(size = 0.75, na.rm = TRUE, alpha = 0.3) +
    # Shade area of mean +/- sd
    ggplot2::stat_summary(
      ggplot2::aes(x = wd_bin), fun.data = "mean_se", geom = "ribbon",
      na.rm = TRUE, alpha = 0.5, fill = "steelblue"
    ) +
    # Add means as points connected by lines
    ggplot2::stat_summary(
      ggplot2::aes(x = wd_bin), fun.y = "mean", geom = "point",
      na.rm = TRUE, color = "orangered"
    ) +
    ggplot2::stat_summary(
      ggplot2::aes(x = wd_bin), fun.y = "mean", geom = "line",
      na.rm = TRUE, color = "orangered"
    ) +
    # Add overall mean as line
    ggplot2::geom_hline(
      data = data %>% dplyr::summarize(mean = mean(!!x, na.rm = TRUE)),
      ggplot2::aes(yintercept = mean), color = "orangered", linetype = 2
    ) +
    # Label x-axis with degrees
    ggplot2::scale_x_continuous(
      breaks = c(seq(0, 360, by = 100)),
      labels = c("0°", "100°", "200°", "300°")
    ) +
    ggplot2::labs(x = NULL, y = paste0(name, " (", units, ")")) +
    ggplot2::theme_bw()
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param timestamp
#' @param step
#' @param remove_outliers
#' @param z
#' @param show_light
#' @param light
#' @param day_thr
#' @param show_summary
#' @param groups
#'
#' @return
#' @export
#'
#' @examples
plot_diurnal <- function(data, x, timestamp = timestamp,
                         step = c("hour", "halfhour"), remove_outliers = TRUE,
                         z = 7, show_light = TRUE, light = NULL, day_thr = 12,
                         show_summary = FALSE, groups = c("none", "month")) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  timestamp <- rlang::enquo(timestamp)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()
  if (!is.null(light)) {
    light <- rlang::enquo(light)
    light_units <- data %>% dplyr::pull(!!light) %>% attr("units")
    light_name <- data %>% dplyr::select(!!light) %>% names()
  }
  groups <- match.arg(groups)
  step <- match.arg(step)
  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n())
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n() - n_all) %>%
      unlist()
  }

  # Pre-condition data for plotting
  data <- data %>%
    # Extract time components from timestamp
    dplyr::mutate(
      hour = lubridate::hour(!!timestamp),
      month = lubridate::month(!!timestamp, label = TRUE),
      month = factor(month)
    ) %>%
    # Specify half-hours if indicated
    {if (step == "halfhour") {
      dplyr::mutate(., hour = hour + lubridate::minute(!!timestamp) / 60)
    } else .} %>%
    # Calculate y-axis scaling factor if light is included
    {if (show_light & groups == "none") {
      dplyr::mutate(., light = normalize(!!light, !!x, hour))
    } else .} %>%
    {if (show_light & groups == "month") {
      dplyr::mutate(., light = normalize(!!light, !!x, month, hour))
    } else .} %>%
    # Group by month if indicated
    {if (groups == "month") {
      dplyr::group_by(., month)
    } else .}

  # Build plot caption according to how data were processed
  if (remove_outliers) {
    cap1 <- paste0(n_removed, " outliers removed using MAD with z = ", z)
  } else cap1 <- NULL
  if (show_light) {
    cap2 <- "\ndaylight hours and incoming radiation shown in orange"
  } else cap2 <- NULL
  if (any(!is.null(cap1), !is.null(cap2))) {
    caption <- paste0(c(cap1, cap2), collapse = "; ")
  } else caption <- NULL

  # Create diurnal averages plot
  out <- data %>%
    # Hack fix to faceting issue: set NA to 0 if a month contains entirely NA
    {if (groups == "month") {
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, dplyr::n_distinct(!!x) < 3, 0)
      )
    } else .} %>%
    #dplyr::filter(!is.na(!!x)) %>%
    ggplot2::ggplot(ggplot2::aes(hour, !!x)) +
    ggplot2::stat_summary(
      fun.data = "mean_se", geom = "ribbon", na.rm = TRUE, alpha = 0.3
    ) +
    ggplot2::stat_summary(fun.y = "mean", geom = "line", na.rm = TRUE) +
    {if (show_light) {
      ggplot2::stat_summary(
        data = data,
        ggplot2::aes(y = light), fun.y = "mean",
        geom = "line", na.rm = TRUE, color = "orange", linetype = 2
      )
    }} +
    ggplot2::scale_x_continuous(
      breaks = c(seq(0, 18, by = 6)),
      labels = c("0:00", "6:00", "12:00", "18:00")
    ) +
    {if (show_light) {
      ggplot2::geom_rect(
        data = data %>%
          {if (groups == "month") {
            dplyr::group_by(., month)
          } else .} %>%
          dplyr::summarize(
            sunrise = mean_suntimes(hour, !!light, day_thr)[1],
            sunset = mean_suntimes(hour, !!light, day_thr)[2]
          ),
        ggplot2::aes(xmin = sunrise, xmax = sunset, ymin = -Inf, ymax = Inf),
        fill = "orange", alpha = 0.2, inherit.aes = FALSE, na.rm = TRUE
      )
    }} +
    {if (show_light) {
      ggplot2::geom_vline(
        data = data %>%
          {if (groups == "month") {
            dplyr::group_by(., month)
          } else .} %>%
          dplyr::summarize(sunrise = mean_suntimes(hour, !!light, day_thr)[1]),
        ggplot2::aes(xintercept = sunrise), color = "orange", na.rm = TRUE
      )
    }} +
    {if (show_light) {
      ggplot2::geom_vline(
        data = data %>%
          {if (groups == "month") {
            dplyr::group_by(., month)
          } else .} %>%
          dplyr::summarize(sunset = mean_suntimes(hour, !!light, day_thr)[2]),
        ggplot2::aes(xintercept = sunset), color = "orange", na.rm = TRUE
      )
    }} +
    ggplot2::labs(
      x = NULL, y = paste0(name, " (", units, ")"), caption = caption
    ) +
    # Facet if grouped by month
    {if (groups == "month") {
      ggplot2::facet_wrap(~month)
    }} +
    ggplot2::theme_bw()
  out
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param timestamp
#' @param remove_outliers
#' @param z
#' @param multiyear
#' @param smooth
#' @param n
#' @param ...
#'
#' @return
#' @export
#'
#' @examples
plot_heatmap <- function(data, x, timestamp = timestamp, remove_outliers = TRUE,
                         z = 7, multiyear = FALSE, smooth = TRUE, n = 5, ...) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  timestamp <- rlang::enquo(timestamp)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()
  # Remove outliers to improve color scale (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all <- data %>% dplyr::filter(is.na(!!x)) %>% dplyr::summarize(n())
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(n() - n_all) %>%
      unlist()
  }
  # Pre-process data for plotting
  data <- data %>%
    # Smooth using moving average (if requested)
    {if (smooth) {
      dplyr::mutate(
        ., !!rlang::as_name(x) := zoo::rollapply(
          !!x, n, mean, na.rm = TRUE, fill = NA
          )
      )
    } else .} %>%
    # Extract time components from timestamp
    dplyr::mutate(
      # Internally shift timestamp -15 minutes to keep all in same year
      !!rlang::as_name(timestamp) := !!timestamp - 900,
      date = lubridate::date(!!timestamp),
      month = lubridate::month(!!timestamp, label = TRUE),
      day = lubridate::day(!!timestamp),
      hour = lubridate::hour(!!timestamp),
      year = lubridate::year(!!timestamp)
    )
  # Build plot caption according to how data were processed
  if (remove_outliers) {
    cap1 <- paste0(n_removed, " outliers removed using MAD with z = ", z)
  } else cap1 <- NULL
  if (smooth) {
    cap2 <- paste0("data smoothed using moving average with n = ", n)
  } else cap2 <- NULL
  if (any(!is.null(cap1), !is.null(cap2))) {
    caption <- paste0(c(cap1, cap2), collapse = "; ")
  } else caption <- NULL
  # Create heatmap plot
  out <- data %>%
    ggplot2::ggplot(ggplot2::aes(day, hour, fill = !!x)) +
    #ggplot2::geom_tile(size = 0.0001, color = "white", ...) +
    ggplot2::geom_raster(interpolate = TRUE, alpha = 0.8, ...) +
    {if (multiyear) {
      ggplot2::facet_grid(year~month)
    } else ggplot2::facet_grid(~month)} +
    ggplot2::labs(caption = caption) +
    ggplot2::labs(x = NULL, y = NULL) +
    ggplot2::scale_fill_distiller(
      palette = "Spectral", na.value = "grey95",
      breaks = function(x) {c(x[1] + .15 * abs(x[2]), x[2] - .2 * abs(x[2]))},
      labels = c("-", "+"),
      guide = ggplot2::guide_colorbar(
        direction = "horizontal", ticks = FALSE, barwidth = 3, barheight = 0.8,
        title = NULL, label.vjust = 5, frame.colour = "white", bins = 30,
        label.theme = element_text(size = 10, face = "bold"),
        frame.linewidth = 2
      )
    ) +
    ggplot2::scale_y_continuous(
      trans = "reverse", breaks = c(seq(0, 18, by = 6)),
      labels = c("0:00", "6:00", "12:00", "18:00"), expand = c(0.005, 0.005)
    ) +
    ggplot2::scale_x_continuous(breaks = c(7, 15, 24), expand = c(0, 0)) +
    ggplot2::theme_bw() +
    ggplot2::theme(
      axis.text.x = ggplot2::element_blank(),
      axis.ticks.x = ggplot2::element_blank(),
      #legend.position = "none",
      legend.key = ggplot2::element_rect(
        fill = "transparent", colour = "transparent"
      ),
      legend.background = ggplot2::element_rect(
        fill = "transparent", colour = "transparent"
      ),
      legend.position = c(0.99, 0),
      legend.justification = c(0.99, 0),
      legend.margin = ggplot2::margin(0, 0, 0, 0),
      legend.box.margin = ggplot2::margin(0, 0, 0, 0),
      legend.title = ggplot2::element_blank(),
      legend.spacing.x = unit(0.1, "cm"),
      panel.border = ggplot2::element_blank(),
      panel.spacing.x = unit(2, "pt"),
      strip.text.x = ggplot2::element_text(
        margin = ggplot2::margin(2, 0, 2, 0)
      ),
      strip.background.x = ggplot2::element_rect(fill = "grey90")
    )
  out
}

## =============================================================================
#' Title
#'
#' @param data
#' @param x
#' @param timestamp
#' @param means
#' @param sd
#' @param remove_outliers
#' @param z
#' @param geom
#' @param smooth
#' @param n
#' @param sample
#' @param perc
#'
#' @return
#' @export
#'
#' @examples
plot_timeseries <- function(data, x, timestamp = timestamp,
                            means = c("none", "daily", "weekly", "monthly"),
                            sd = TRUE, remove_outliers = TRUE, z = 7,
                            geom = c("point", "line"), smooth = TRUE, n = 5,
                            sample = FALSE, perc = 0.5) {
  # Get variables and attributes
  x <- rlang::enquo(x)
  timestamp <- rlang::enquo(timestamp)
  units <- data %>% dplyr::pull(!!x) %>% attr("units")
  name <- data %>% dplyr::select(!!x) %>% names()
  geom <- match.arg(geom)
  means <- match.arg(means)

  # Extract time components from timestamp
  data <- data %>%
    dplyr::mutate(
      year = lubridate::year(!!timestamp),
      date = lubridate::date(!!timestamp),
      week = lubridate::week(!!timestamp) * 7,
      week = paste0(year, "-", week),
      month = format(!!timestamp, "%Y-%m")
    )

  # Remove outliers to improve plot range (if requested)
  if (remove_outliers) {
    med <- data %>%
      dplyr::summarize(median(!!x, na.rm = TRUE)) %>% unlist() %>% unname()
    n_all <- data %>% dplyr::filter(is.na(!!x)) %>% dplyr::summarize(n())
    data <- data %>%
      dplyr::group_by(!!x > med) %>%
      dplyr::mutate(
        ., !!rlang::as_name(x) := replace(!!x, get_outliers(!!x, z = z), NA)
      ) %>%
      dplyr::ungroup()
    # Calculate number of outliers removed
    n_removed <- data %>%
      dplyr::filter(is.na(!!x)) %>%
      dplyr::summarize(dplyr::n() - n_all) %>%
      unlist()
  }

  # Build plot caption according to how data were processed
  if (remove_outliers) {
    cap1 <- paste0(n_removed, " outliers removed using MAD with z = ", z)
  } else cap1 <- NULL
  if (smooth) {
    cap2 <- paste0("means smoothed using moving average with n = ", n)
  } else cap2 <- NULL
  if (any(!is.null(cap1), !is.null(cap2))) {
    caption <- paste0(c(cap1, cap2), collapse = "; ")
  } else caption <- NULL

  # Create time series plot
  out <- data %>%
    # Sample to avoid overplotting if specified
    {if (sample) {
      dplyr::sample_frac(., perc)
    } else .} %>%
    ggplot2::ggplot(ggplot2::aes(!!timestamp, !!x)) +
    {if (geom == "point") {
      ggplot2::geom_point(alpha = 0.3, size = 1, stroke = 0.2, na.rm = TRUE)
    }} +
    {if (geom == "line") {
      ggplot2::geom_line(alpha = 0.8, size = 0.5, na.rm = TRUE)
    }} +
    ggplot2::labs(
      x = NULL, y = paste0(name, " (", units, ")"), caption = caption
    ) +
    ggplot2::theme_bw()

  # Add aggregated means if specified
  if (means != "none") {
    grouped <- data %>%
    {if (means == "daily") {
      dplyr::mutate(., group = date)
    } else .} %>%
    {if (means == "weekly") {
      dplyr::mutate(., group = week)
    } else .} %>%
    {if (means == "monthly") {
      dplyr::mutate(., group = month)
    } else .} %>%
      dplyr::group_by(group) %>%
      dplyr::summarize(
        nrow = sum(!is.na(!!x)),
        mean = dplyr::if_else(nrow > 5, mean(!!x, na.rm = TRUE), NA_real_),
        sd = dplyr::if_else(nrow > 5, sd(!!x, na.rm = TRUE), NA_real_)
      ) %>%
      {if (smooth) {
        dplyr::mutate(
          ., mean = zoo::rollapply(mean, n, mean, na.rm = TRUE, fill = NA),
          sd = zoo::rollapply(sd, n, mean, na.rm = TRUE, fill = NA)
        )
      } else .} %>%
      {if (means == "daily") {
        dplyr::mutate(
          ., group = lubridate::parse_date_time(group, "Ymd") + 43200
        )
      } else .} %>%
      {if (means == "weekly") {
        dplyr::mutate(
          ., group = lubridate::parse_date_time(group, "Yj") + 302400
        )
      } else .} %>%
      {if (means == "monthly") {
        dplyr::mutate(
          ., group = lubridate::parse_date_time(group, "Ym") + 1252800)
      } else .}

    # Add sd components to plot if specified
    if (sd) {
      out <- out +
        # Shade area of mean +/- sd
        ggplot2::geom_ribbon(
          data = grouped,
          ggplot2::aes(x = group, ymin = mean - sd, ymax = mean + sd),
          fill = "steelblue", alpha = 0.5, inherit.aes = FALSE
        ) +
        # Bottom shading border
        ggplot2::geom_line(
          data = grouped, ggplot2::aes(group, mean - sd),
          color = "steelblue", na.rm = TRUE
        ) +
        # Top shading border
        ggplot2::geom_line(
          data = grouped, ggplot2::aes(group, mean + sd),
          color = "steelblue", na.rm = TRUE
        )
    }
    # Add means as lines
    out <- out +
      ggplot2::geom_line(
        data = grouped, ggplot2::aes(group, mean),
        color = "orangered", na.rm = TRUE, size = 0.75
      )
  }
  out
}


## =============================================================================
#' Plot Precheck Time Series of Eddy Covariance Data
#'
#' @param data
#' @param var
#' @param ylim_list
#'
#' @return A half-hourly time series plot.
plot_precheck <- function(data, var, ylim_list = NULL) {
  if (!is.null(ylim_list)) {
    # Define custom data range if provided
    ylim <- ylim_list[[var]]
  } else {
    # Define automatic data range to exclude outliers
    med <- median(data[, var], na.rm = TRUE)
    iqr <- IQR(data[, var], na.rm = TRUE)
    ylim <- c(med - 10 * iqr, med + 10 * iqr)
  }
  range <- range(data[, var], na.rm = TRUE) # original data range
  # Return plot range to original range if there are no outliers
  if (ylim[1] < range[1]) ylim[1] <- range[1]
  if (ylim[2] > range[2]) ylim[2] <- range[2]
  omit <- data[data[, var] < ylim[1] | data[, var] > ylim[2], var]
  omit <- omit[!is.na(omit)] # data outliers
  range_omit <- range(omit) # range of outliers
  min_omit <- range(omit)[1] # min outlier value
  max_omit <- range(omit)[2] # max outlier value
  n_omit <- length(omit) # number of outliers
  if (any(!is.finite(range))) {
    min_omit <- range[1]
    max_omit <- range[2]
  } else {
    if (min_omit > range[1]) min_omit <- range[1]
    if (max_omit < range[2]) max_omit <- range[2]
  }
  plot(
    data$timestamp, data[, var],
    ylim = ylim,
    pch = 20, cex = 0.5, col = adjustcolor("black", 0.5),
    xaxt = "n", xlab = "timestamp", ylab = var, main = var
  )
  r <- as.POSIXct(round(range(data$timestamp), "months"))
  axis.POSIXct(1, at = seq(r[1], r[2], by = "month"), format = "%b-%y")
  mtext(paste(
    n_omit, "outliers masked; original range", signif(min_omit, 5),
    "to", signif(max_omit, 5)
  ), side = 3)
}

# ==============================================================================
#' Set Range for Plotting
#'
#' Ensure that the range extracted from \code{x} is valid for plotting.
#'
#' Set range for a given numeric vector \code{x} subsetted by a logical vector
#' \code{filter}. If the subset contains any finite value, finite range of the
#' \code{x} subset is returned. Else if \code{x} contains any finite value,
#' finite range of \code{x} is returned. When no finite value can be found in
#' \code{x}, range is set manually (default \code{man = c(0, 0)}).
#'
#' If \code{x} length is higher than \code{filter} length, \code{filter} will be
#' recycled.
#'
#' @param x A numeric vector.
#' @param filter A logical vector that can be recycled if necessary to match the
#'   length of \code{x}.
#' @param man A numeric vector of length 2.
#'
#' @seealso \code{\link{range}}.
#'
#' @examples
#' (aa <- c(1, NA, Inf, -Inf, 3, NaN, rep(NA, 4)))
#' range(aa, finite = TRUE)
#' setRange(aa, TRUE) # same effect
#' # Useful when applying filters
#' aa[rep(c(F, T), each = 5)]
#' suppressWarnings(range(aa[rep(c(F, T), each = 5)], finite = TRUE))
#' setRange(aa, rep(c(F, T), each = 5)) # range taken from unfiltered 'aa'
#' setRange(aa[c(FALSE, TRUE)]) # No finite values in 'x', applies 'man' range
set_range <- function(x = NA, filter = TRUE, man = c(0, 0), mad_adj = FALSE) {
  if (!is.numeric(man) || length(man) != 2 || any(!is.finite(man))) {
    stop("'man' must be numeric vector with 2 finite values.", call. = FALSE)
  }
  if (is.null(x) || !is.numeric(x) && !all(is.na(x))) {
    stop("'x' must be a vector containing numeric or NA values.", call. = FALSE)
  }
  if (!is.logical(filter)) {
    stop("'filter' must be a logical vector.", call. = FALSE)
  }
  if (length(x) %% length(filter) != 0) {
    warning(paste(
      "'x' length [", length(x), "] is not a multiple of 'filter'",
      " length [", length(filter), "].",
      sep = ""
    ), call. = FALSE)
  }
  if (any(is.finite(x[filter]))) {
    x <- x[filter]
    if (mad_adj == TRUE) {
      med <- median(x, na.rm = TRUE)
      mad <- mad(x, na.rm = TRUE)
      x[x > med + 7 * mad | x < med - 7 * mad] <- NA
    }
    return(range(x, finite = TRUE))
  } else if (any(is.finite(x))) {
    if (mad_adj == TRUE) {
      med <- median(x, na.rm = TRUE)
      mad <- mad(x, na.rm = TRUE)
      x[x > med + 7 * mad | x < med - 7 * mad] <- NA
    }
    return(range(x, finite = TRUE))
  } else {
    return(man)
  }
}

# ==============================================================================
#' Plot Time Series of Eddy Covariance Data
#'
#' Visualize flux measurements together with micrometeorological or other
#' related variables in monthly and weekly intervals. Missing values, scales of
#' axes and plotting regions are treated in an automated way.
#'
#' The data frame \code{x} is expected to have certain properties. It is
#' required that it contains column named \code{"timestamp"} of class
#' \code{"POSIXt"} with regular sequence of date-time values, typically with
#' (half-)hourly frequency. Any missing values in \code{"timestamp"} are not
#' allowed. Thus, if no records exist for given date-time value, it still has to
#' be included. It also has to contain required (depends on the argument values
#' and applied modules) column names. If required auxiliary variable is not
#' available (e.g. \code{"P"}, \code{"T"}), it still has to be included in
#' \code{x} as a named column with \code{NA} values. The \code{x} column defined
#' by argument \code{flux} is the only one that has to contain at least one
#' non-missing value.
#'
#' If \code{skip = "weekly"}, minimum requirements for \code{x} column names are
#' \code{"timestamp"} and \code{flux}. If \code{skip = "none"} or
#' \code{"monthly"}, \code{"P"} and respective names specified in 'Modules' (see
#' below) are also required.
#'
#' Variable names for plot labels are extracted from required column names of
#' \code{x}. Units are extracted from \code{x} if they are present as
#' \code{units} attributes of required columns. If missing, \code{"-"} string is
#' used instead.
#'
#' Plotting is separated into two stages. Firstly, \code{flux} time-series data
#' are drawn in monthly intervals. Monthly plotting regions consist of four
#' figures on top of each other representing separately four consecutive months.
#' Secondly, if \code{skip = "none"} or \code{"monthly"}, \code{flux}
#' time-series data are drawn together with auxiliary data in weekly intervals.
#' Weekly plotting regions are described in 'Modules' section (see below).
#'
#' @section Modules: Applies only if \code{skip = "none"} or \code{"monthly"}.
#'   Plotting in weekly intervals is simplified by using predefined modules.
#'   Their main purpose is to achieve well-arranged display of auxiliary
#'   variables. Weekly plotting regions consist of two figures representing
#'   separately two consecutive weeks. Each figure contains three panels on top
#'   of each other. The middle panel always contains the values from \code{flux}
#'   and \code{"P"} columns of \code{x}. Variables used in the upper and lower
#'   panel can be changed by \code{panel_top} and \code{panel_bottom}. These
#'   arguments specify the respective modules that will be loaded (can be the
#'   same) and thus also a certain set of required column names of \code{x}
#'   (variables).
#'
#'   Available modules are: \itemize{ \item T_light: requires \code{"Tair"},
#'   \code{"Tsoil"} and selected \code{light} columns. \item VPD_Rn: requires
#'   \code{"VPD"} and \code{"Rn"} columns. \item H_err_var: requires
#'   \code{"rand_err_H"} and \code{"ts_var"} columns.}
#'
#' @section Quality Control: \code{qc_flag} and \code{test} relate to QC flags
#'   available for the specified \code{flux}. Only QC scheme using QC flag range
#'   0 - 2 is supported.
#'
#'   Flags specified by \code{qc_flag} separate corresponding flux values to two
#'   groups. \emph{Used data} for flags 0 and 1 (black points) and
#'   \emph{Excluded data} for flags 2 and \code{NA} (grey points). The range of
#'   y-axis limit in the figures with \code{flux} values is based on \emph{Used
#'   data} only. If \code{qc_flag = "none"}, all data are \emph{Used data}. The
#'   time interval used for setting the range is one month both for monthly
#'   (respective month) and weekly (month centered on the respective week)
#'   plots.
#'
#'   In order to emphasize \code{flux} values with lower quality, \code{test}
#'   can be specified. Values with QC flag = 1 have green center of the point.
#'   Values with QC flag = 2 or \code{NA} have red center of the point.
#'
#'   NB: \code{flux} data with \code{NA} values are always \emph{Excluded data}
#'   and cannot be emphasized (\code{NA} values are not drawn).
#'
#' @section Gap-filling and NEE separation: Gap-filled flux values can be
#'   displayed using \code{flux_gf} as a line overlaying the \code{flux} values.
#'   If \code{NEE_sep = TRUE}, columns \code{"Reco"} and \code{"GPP"} are
#'   expected in \code{x}. \code{GPP_check} is evaluated only if \code{NEE_sep =
#'   TRUE}. If \code{GPP_check = TRUE}, mean GPP is checked to be negative so it
#'   compares well with NEE and corrected accordingly if needed. Thus negative
#'   values denote carbon sink. Such test can fail for measurements above
#'   surfaces without vegetation or ecosystems with insignificant sink. In that
#'   case set \code{GPP_check = FALSE} and check and correct GPP sign convention
#'   manually.
#'
#' @section Abbreviations: \itemize{ \item H: Sensible heat flux [W m-2] \item
#'   NEE: Net Ecosystem Exchange [umol m-2 s-1] \item GPP: Gross Primary
#'   Production [umol m-2 s-1] \item Reco: Ecosystem Respiration [umol m-2 s-1]
#'   \item QC: Quality Control \item P: Precipitation [mm] \item PAR:
#'   Photosynthetic Active Radiation [umol m-2 s-1] \item GR: Global Radiation
#'   [W m-2] \item T: Temperature [degC] \item Tair: Air Temperature [degC]
#'   \item Tsoil: Soil Temperature [degC] \item VPD: Vapor Pressure Deficit
#'   [hPa] \item Rn: Net Radiation [W m-2] \item rand_err_H: random error of H
#'   [W m-2]; in plots abbreviated as H_re \item ts_var: sonic temperature
#'   variance [K2]}
#'
#' @param x A data frame with column names representing required variables. See
#'   'Details' below.
#' @param flux A character string. Specifies the column name in \code{x} with
#'   flux values.
#' @param qc_flag A character string. Specifies the column name in \code{x} with
#'   flux related quality control flag used for visualisation of data quality
#'   and setting of y-axis range. If "none" is provided, all data will be used.
#'   See 'Quality Control'.
#' @param test A character string. Specifies the column name in \code{x} with
#'   quality control flag for visualisation of its effect on the data. If "none"
#'   is provided, no visualisation will be done. See 'Quality Control'.
#' @param flux_gf A character string. Specifies the column name in \code{x} with
#'   gap-filled flux values.
#' @param NEE_sep A logical value. Determines whether NEE separation should be
#'   visualized. If \code{TRUE}, columns \code{"Reco"} and \code{"GPP"} are
#'   expected in \code{x}.
#' @param skip A character string. Determines whether plotting should be done in
#'   monthly (\code{skip = "weekly"}), weekly intervals (\code{skip =
#'   "monthly"}) or both (\code{skip = "none"}).
#' @param panel_top A character string. Selects one of the available modules for
#'   plotting additional variables. This module is displayed above the panel
#'   with fluxes in weekly plots. Can be abbreviated. See 'Modules'.
#' @param panel_bottom A character string. Selects one of the available modules
#'   for plotting additional variables. This module is displayed below the panel
#'   with fluxes in weekly plots. Can be abbreviated. See 'Modules'.
#' @param light A character string. Required only for the \code{"T_light"}
#'   module. Selects preferred variable for incoming light intensity.
#'   \code{"PAR"} or \code{"GR"} is allowed. Can be abbreviated.
#' @param GPP_check A logical value. If \code{TRUE}, column \code{"GPP"} is
#'   checked to be of the same sign convention as \code{"NEE"} column and
#'   corrected in the case it is not. Required only if \code{NEE_sep = TRUE}.
#' @param document A logical value. If \code{TRUE}, values of \code{qc_flag} and
#'   \code{test} arguments are documented in both monthly and weekly plots.
#'
#' @seealso \code{\link{read_eddy}} and \code{\link{strptime_eddy}}.
plot_eddy <- function(x, flux, qc_flag = "none", test = "none",
                      flux_gf = "none", NEE_sep = FALSE,
                      skip = c("none", "monthly", "weekly"),
                      panel_top = c("t_light", "vpd_rn", "h_err_var"),
                      panel_bottom = c("vpd_rn", "t_light", "h_err_var"),
                      light = c("ppfd", "rg"), GPP_check = TRUE,
                      document = TRUE, mad_adj_lims = FALSE) {
  x_names <- names(x)
  check_input(x, "data_frame")
  if (any(!sapply(list(flux, qc_flag, test), is.character))) {
    stop("'flux', 'qc_flag', 'test' must be of class character.")
  }
  req_vars <- c("timestamp", flux)
  req <- !c(qc_flag, test, flux_gf) %in% "none"
  req_vars <- c(req_vars, c(qc_flag, test, flux_gf)[req])
  if (NEE_sep) req_vars <- c(req_vars, "Reco", "GPP")
  skip <- match.arg(skip)
  if (skip != "weekly") {
    req_vars <- c(req_vars, "p")
    panel_top <- match.arg(panel_top)
    panel_bottom <- match.arg(panel_bottom)
    if ("t_light" %in% c(panel_top, panel_bottom)) {
      light <- match.arg(light)
      req_vars <- c(req_vars, light, "tair", "tsoil")
    }
    if ("vpd_rn" %in% c(panel_top, panel_bottom)) {
      req_vars <- c(req_vars, "vpd", "rn")
    }
    if ("h_err_var" %in% c(panel_top, panel_bottom)) {
      req_vars <- c(req_vars, "unc_h", "ts_sd")
    }
  }
  if (!all(req_vars %in% x_names)) {
    stop(paste("Missing", paste0(
      req_vars[!(req_vars %in% x_names)], collapse = ", "
    ), "."))
  }
  if (!inherits(x$timestamp, "POSIXt")) {
    stop("'x$timestamp' must be of class 'POSIXt'.")
  }
  num_vars <- req_vars[!req_vars %in% c("timestamp", qc_flag, test)]
  check1 <- function(x) {
    !is.numeric(x) && !all(is.na(x))
  }
  res_c1 <- sapply(x[, num_vars], check1)
  if (any(res_c1)) {
    stop(paste0(
      "Columns [", paste0(num_vars[res_c1], collapse = ", "),
      "] in 'x' must be of class numeric or contain NAs only."
    ))
  }
  if (qc_flag != "none" || test != "none") {
    check2 <- function(x) {
      !is.numeric(x) && !is.logical(x)
    }
    res_c2 <- sapply(x[, c(qc_flag, test)[c(qc_flag, test) != "none"]], check2)
    if (any(res_c2)) {
      stop(paste0(
        "Columns [", paste0(c(qc_flag, test)[res_c2], collapse = ", "),
        "] in 'x' must be of class numeric or logical."
      ))
    }
  }
  time <- as.POSIXlt(x$timestamp)
  if (any(is.na(time))) stop("NAs in 'timestamp' not allowed.")
  if (any(diff(as.numeric(time)) != mean(diff(as.numeric(time))))) {
    stop("Timestamp does not form regular sequence.")
  }
  units <- openeddy::units(x, names = TRUE)
  wrap <- function(x) paste("[", x, "]", sep = "")
  date <- as.Date(x$timestamp)
  vals <- x[, flux]
  qc <- if (qc_flag == "none") 0L else x[, qc_flag] # qc_flag 0: show all data
  qc[is.na(qc)] <- 2L # NA qc is interpreted as flag 2
  exalt <- if (test == "none") 0L else x[, test] # exalt = 0: exalt not used
  exalt[is.na(exalt)] <- 2L # NA qc is interpreted as flag 2
  if (flux_gf != "none") vals_gf <- x[, flux_gf]
  if (NEE_sep) {
    Reco <- x[, "Reco"]
    GPP <- x[, "GPP"]
    if (GPP_check) if (mean(GPP, na.rm = TRUE) > 0) GPP <- -GPP
  }
  use <- qc < 2 & !is.na(vals)
  if (!any(use)) {
    stop("No non-NA values with accepted quality in 'flux'.")
  }
  # Correcting $year (counting since 1900) and $mon (0-11 range)
  # Selects first day of month in the dataset
  day1 <- as.Date(paste0(time$year[1] + 1900, "-", time$mon[1] + 1, "-01"))
  # Graphical display of data in monthly periods ===============================
  if (skip != "monthly") {
    # grey = excluded data, green = used data
    # Number of intervals with right side closure (+1)
    n_int <- length(unique(paste(time$year, time$mon))) + 1L
    # Create monthly intervals for plotting
    int <- seq(from = day1, length.out = n_int, by = "1 month")
    op <- par(mfcol = c(4, 1), mar = c(3, 0, 0, 0), oma = c(2, 6, 1, 1))
    for (i in 1:(n_int - 1L)) {
      mon <- date >= int[i] & date < int[i + 1L]
      # keep the lenght of time[mon] but remove excluded data
      # (needed for lines)
      showVals <- vals[mon]
      showVals[!use[mon]] <- NA
      # xaxis has to be one day longer to simplify plotting
      xaxis <- seq(from = int[i], to = int[i + 1L], by = "1 day")
      xaxis <- as.POSIXct(strptime(xaxis, "%Y-%m-%d", tz = "GMT"))
      y <- vals[use]
      f <- mon[use]
      if (flux_gf != "none") {
        y <- c(y, vals_gf)
        f <- c(f, mon)
      }
      if (NEE_sep) {
        y <- c(y, Reco, GPP)
        f <- c(f, rep(mon, 2))
      }
      plot(
        time[mon], vals[mon], type = "n", xaxt = "n", yaxt = "n",
        ylab = "", xlab = "", panel.first = grid(nx = NA, ny = NULL),
        xlim = range(xaxis), ylim = set_range(y, f, mad_adj = mad_adj_lims)
      )
      # Halfday shift of ticks to mark the middle of day
      axis.POSIXct(
        1, at = seq(min(xaxis), max(xaxis), by = "5 days") + 43200,
        format = "%Y/%m/%d", padj = -0.5
      )
      axis(2)
      abline(v = xaxis, col = "grey", lty = 3)
      lines(
        time[mon], vals[mon], type = "o", pch = 19, cex = 0.75, col = "grey"
      )
      lines(time[mon], showVals, type = "o", pch = 19, cex = 0.75)
      if (flux_gf != "none") {
        lines(
          time[mon], vals_gf[mon],
          type = "l", pch = 19, cex = 0.4, col = "forestgreen"
        )
      }
      if (NEE_sep) {
        lines(time[mon], Reco[mon], col = "firebrick3")
        lines(time[mon], GPP[mon], col = "dodgerblue3")
      }
      points(
        time[mon & (exalt == 1)], vals[mon & (exalt == 1)],
        col = "green", pch = 19, cex = 0.4
      )
      points(
        time[mon & (exalt == 2)], vals[mon & (exalt == 2)],
        col = "red", pch = 19, cex = 0.4
      )
      mtext(wrap(units[flux]), 2, line = 2.2, cex = 0.8)
      mtext(flux, 2, line = 3.4, cex = 1.2)
      if (i %% 4 == 1) {
        legend(
          "topleft", pch = 19, lty = c(1, 1, NA, NA), bty = "n",
          legend = c(
            "Used data", "Excluded data", "Test flag = 1", "Test flag = 2"
          ),
          col = c("black", "grey", "green", "red")
        )
        if (flux_gf != "none" || NEE_sep) {
          legend(
            "topright", lwd = 2, bty = "n",
            legend = c("Gap-filled data", if (NEE_sep) c("Reco", "GPP")),
            col = c(
              "forestgreen", if (NEE_sep) c("firebrick3", "dodgerblue3")
            )
          )
        }
      }
      if (document && i %% 4 == 3) {
        mtext(
          paste0("qc_flag = '", qc_flag, "', test = '", test, "'"),
          side = 3, line = 0, adj = 0, cex = 0.6)
      }
      if (i %% 4 == 0) mtext("Day of year", 1, line = 3, cex = 1.2)
    }
    par(op)
  }
  if (skip != "weekly") {
    p <- x$p
    # Modules ==================================================================
    modules <- list()
    if ("t_light" %in% c(panel_top, panel_bottom)) {
      ppfd <- x[, light]
      tair <- x$tair
      tsoil <- x$tsoil
      modules$t_light <- function() {
        plot(
          time[week], tair[week], xlim = range(xaxis),
          ylim = set_range(c(tair, tsoil), mon), type = "n",
          panel.first = grid(nx = NA, ny = NULL), xaxt = "n", yaxt = "n",
          ylab = "", xlab = ""
        )
        abline(v = xaxis, col = "grey", lty = 3)
        axis(2, line = 1.8, padj = 0.9, tcl = -0.3)
        lines(time[week], tair[week], lwd = 2, col = "dodgerblue")
        lines(time[week], tsoil[week], lwd = 2, col = "red1")
        par(new = TRUE)
        plot(
          time[week], ppfd[week], xlim = range(xaxis),
          ylim = set_range(ppfd, mon), type = "l", col = "gold", lwd = 2,
          xaxt = "n", yaxt = "n", xlab = "", ylab = ""
        )
        axis(4, padj = -0.9, tcl = -0.3)
        mtext(paste("T", wrap(units["tair"])), 2, line = 3.6)
        mtext(light, 4, line = 2)
        mtext(wrap(units[light]), 4, line = 3.6, cex = 0.8)
        if (i %% 2 == 1) {
          legend(
            "topleft", legend = c("tair", "tsoil", light), lty = 1, lwd = 2,
            col = c("dodgerblue", "red1", "gold"), bty = "n"
          )
        }
      }
    }
    if ("vpd_rn" %in% c(panel_top, panel_bottom)) {
      vpd <- x$vpd
      rn <- x$rn
      modules$vpd_rn <- function() {
        plot(
          time[week], vpd[week], xlim = range(xaxis),
          ylim = set_range(vpd, mon), panel.first = grid(nx = NA, ny = NULL),
          type = "n", xaxt = "n", yaxt = "n", ylab = "", xlab = ""
        )
        abline(v = xaxis, col = "grey", lty = 3)
        axis(2, line = 1.8, padj = 0.9, tcl = -0.3)
        lines(time[week], vpd[week], lwd = 2, col = "mediumorchid")
        par(new = TRUE)
        plot(
          time[week], rn[week], xlim = range(xaxis), ylim = set_range(rn, mon),
          type = "l", col = "chocolate1", lwd = 2,
          xaxt = "n", yaxt = "n", xlab = "", ylab = ""
        )
        axis(4, padj = -0.9, tcl = -0.3)
        mtext(paste("vpd", wrap(units["vpd"])), 2, line = 3.6)
        mtext("rn", 4, line = 2)
        mtext(wrap(units["rn"]), 4, line = 3.6, cex = 0.8)
        if (i %% 2 == 1) {
          legend(
            "topleft", legend = c("vpd", "rn"), lty = 1, lwd = 2,
            col = c("mediumorchid", "chocolate1"), bty = "n"
          )
        }
      }
    }
    if ("h_err_var" %in% c(panel_top, panel_bottom)) {
      unc_h <- x$unc_h
      ts_sd <- x$ts_sd
      modules$H_err_var <- function() {
        plot(
          time[week], unc_h[week],
          xlim = range(xaxis), ylim = set_range(unc_h, mon),
          panel.first = grid(nx = NA, ny = NULL),
          type = "n", xaxt = "n", yaxt = "n", ylab = "", xlab = ""
        )
        abline(v = xaxis, col = "grey", lty = 3)
        axis(2, line = 1.8, padj = 0.9, tcl = -0.3)
        lines(time[week], unc_h[week], lwd = 2, col = "mediumorchid")
        par(new = TRUE)
        plot(
          time[week], ts_sd[week],
          xlim = range(xaxis), ylim = set_range(ts_sd, mon),
          type = "l", col = "chocolate1", lwd = 2,
          xaxt = "n", yaxt = "n", xlab = "", ylab = ""
        )
        axis(4, padj = -0.9, tcl = -0.3)
        mtext(paste("unc_h", wrap(units["unc_h"])), 2, line = 3.6)
        mtext("ts_sd", 4, line = 2)
        mtext(wrap(units["ts_sd"]), 4, line = 3.6, cex = 0.8)
        if (i %% 2 == 1) {
          legend(
            "topleft", legend = c("unc_h", "ts_sd"), lty = 1, lwd = 2,
            col = c("mediumorchid", "chocolate1"), bty = "n"
          )
        }
      }
    }
    # Graphical display of data in weekly periods===============================
    # Number of intervals with right side closure (+1)
    n_int <- length(seq(from = day1, to = date[nrow(x)], by = "1 week")) + 1L
    # Create weekly intervals for plotting
    int <- seq(from = day1, length.out = n_int, by = "1 week")
    # Compute xlim for barplot (60 * 24 = 1440: minutes in day; 7 days in week)
    tdiff <- as.numeric(time[2] - time[1])
    barxaxis <- 1440 / tdiff * 7
    # Only 6 plots are printed, slots 7 and 8 reserved for margins
    panels <- c(1, 1, 1, 2, 2, 2, 3, 3, 3, 7, 7, 4, 4, 4, 5, 5, 5, 6, 6, 6, 8)
    def_par <- par(no.readonly = TRUE)
    par(mar = c(0, 0, 0, 0), oma = c(2.5, 6, 1, 6))
    for (i in 1:(n_int - 1L)) {
      if (i %% 2 == 1) layout(panels)
      week <- date >= int[i] & date < int[i + 1L]
      if (!length(time[week])) next
      center <- int[i] + 3.5
      mon <- date >= (center - 15) & date < (center + 15)
      # xaxis has to be one day longer to simplify plotting
      xaxis <- seq(from = int[i], to = int[i + 1], by = "1 day")
      xaxis <- as.POSIXct(strptime(xaxis, "%Y-%m-%d", tz = "GMT"))
      modules[[panel_top]]()
      if (document && i %% 2 == 0) {
        mtext(
          paste0("qc_flag = '", qc_flag, "', test = '", test, "'"),
          side = 3, line = 0, adj = 0, cex = 0.6
        )
      }
      y <- vals[use]
      f <- mon[use]
      if (flux_gf != "none") {
        y <- c(y, vals_gf)
        f <- c(f, mon)
      }
      if (NEE_sep) {
        y <- c(y, Reco, GPP)
        f <- c(f, rep(mon, 2))
      }
      yRange <- set_range(y, f)
      plot(
        time[week], vals[week], type = "n", xaxt = "n", yaxt = "n",
        ylab = "", xlab = "", panel.first = grid(nx = NA, ny = NULL),
        xlim = range(xaxis), ylim = yRange
      )
      abline(v = xaxis, col = "grey", lty = 3)
      axis(2, padj = 0.9, tcl = -0.3)
      par(new = TRUE)
      barplot(
        p[week], ylim = rev(set_range(p, mon)), xlim = c(0, barxaxis),
        col = "lightskyblue", space = 0, border = NA,
        xaxt = "n", yaxt = "n", xlab = "", ylab = ""
      )
      axis(4, line = 1.8, padj = -0.9, tcl = -0.3)
      mtext(paste("p", wrap(units["p"])), 4, line = 3.6)
      par(new = TRUE)
      plot(
        time[week], vals[week],
        xlim = range(xaxis), ylim = yRange, type = "o", pch = 19, cex = 0.75,
        xaxt = "n", yaxt = "n", xlab = "", ylab = "", col = "grey"
      )
      # Keep length of time[week] but remove excluded data (for plotting lines)
      showVals <- vals[week]
      showVals[!use[week]] <- NA
      lines(
        time[week], showVals, type = "o", col = "black", pch = 19, cex = 0.75
      )
      if (flux_gf != "none") {
        lines(
          time[week], vals_gf[week],
          type = "l", pch = 19, cex = 0.4, col = "forestgreen"
        )
      }
      if (NEE_sep) {
        lines(
          time[week], Reco[week],
          type = "l", pch = 19, cex = 0.4, col = "firebrick3"
        )
        lines(
          time[week], GPP[week],
          type = "l", pch = 19, cex = 0.4, col = "dodgerblue3"
        )
      }
      points(
        time[week & (exalt == 1)], vals[week & (exalt == 1)],
        col = "green", pch = 19, cex = 0.4
      )
      points(
        time[week & (exalt == 2)], vals[week & (exalt == 2)],
        col = "red", pch = 19, cex = 0.4
      )
      mtext(wrap(units[flux]), 2, line = 2, cex = 0.8)
      mtext(flux, 2, line = 3.6)
      if (i %% 2 == 1) {
        legend(
          "topleft", col = c("black", "grey", "green", "red", "lightskyblue"),
          legend = c(
            "Used data", "Excluded data", "Test flag = 1", "Test flag = 2", "P"
          ),
          bty = "n", pch = c(19, 19, 19, 19, 15), lty = c(1, 1, NA, NA, NA)
        )
        if (flux_gf != "none" || NEE_sep) {
          legend(
            "topright",
            legend = c("Gap-filled data", if (NEE_sep) c("Reco", "GPP")),
            col = c("forestgreen", if (NEE_sep) c("firebrick3", "dodgerblue3")),
            lwd = 2, bty = "n"
          )
        }
      }
      modules[[panel_bottom]]()
      # Halfday shift of ticks to mark the middle of day
      axis.POSIXct(
        1, at = seq(min(xaxis), max(xaxis), by = "days") + 12 * 3600,
        format = "%Y/%m/%d", padj = -0.5
      )
      if (i %% 2 == 0) mtext("Day of year", 1, line = 3, cex = 1.2)
    }
    par(def_par)
  }
}

## =============================================================================
#' Plot Daily Sums for a Variable in Specified Year
#'
#' The daily sums for a single year are plotted to the current device, scaled to
#' all data. The daily sums are only calculated for days with complete data.
#' Original name: sEddyProc_sPlotDailySumsY.
#'
#' @param data A data frame with variables to plot.
#' @param var A (gap-filled) variable to plot.
#' @param var_unc The uncertainty estimates for the variable.
#' @param year The year to plot.
#' @param tf Time conversion factor with default per second to per day.
#' @param mf Mass conversion factor with default from umol CO2 to g C.
#' @param unit The unit of the daily sums.
#' @param dts A numeric integer indicating the number of daily time steps (24 or
#'   48.)
#'
#' @details
#' This function first computes the average flux for each day. If the original
#' unit is not "per day", then it need to be converted to "per day" by argument
#' \code{tf}. Furthermore, a change of the mass unit is provided by argument
#' \code{mf}. The default parameters assume original units of umol CO2 / m2 /
#' second and convert to gC / m2 / day. The conversion factors allow plotting
#' variables with different units.
#'
#' @return
#'
#' @export
plot_daily <- function(data, var, var_unc = "none", year, tf = NULL, mf = NULL,
                       unit = "gC/m2/day", dts = 48) {
  # Set plot contents
  timestamp <- data$timestamp
  if (is.null(tf)) tf <- 3600 * 24 else tf <- tf
  if (is.null(mf)) mf <- (44.0096 / 1e6) * (12.011 / 44.0096) else mf <- mf
  data <- apply_qc(data, var, "none", NA)
  full_year <- expand_year(data, timestamp, year, dts)
  data_time <- full_year$timestamp
  data_plot <- full_year$data
  if (var_unc != "none") {
    sd <- apply_qc(data, var_unc, "none", NA)
    data_sd <- expand_year(sd, timestamp, year, dts)$data
  } else {
    # Set uncertainties to zero
    data_sd <- (rep(0, length(data_plot)))
  }
  # Uncertainty only used where data
  data_sd <- ifelse(is.na(data_plot), NA, data_sd)
  # If there is data but no uncertainty estimates, an empty box will be plotted
  missing_unc <- sum(!is.na(data_plot) & is.na(data_sd))
  # Compute daily sums
  daily_recs <- dts
  data_year <- matrix(as.numeric(format(data_time, "%Y")), nrow = dts)[1, ]
  data_doy <- matrix(as.numeric(format(data_time, "%j")), nrow = dts)[1, ]
  data_avg <- (1 / dts) * apply(matrix(data_plot, nrow = daily_recs), 2, mean)
  data_sum <- data_avg * tf * mf
  # Sum of squares function
  sos <- function(x, ...) {
    sum(x^2, ...)
  }
  unc <- apply(matrix(data_sd, nrow = daily_recs), 2, sos)
  data_unc <- (1 / daily_recs) * sqrt(unc) * tf * mf
  # Scale to all data
  y_min <- min(data_sum - data_unc, na.rm = TRUE)
  y_max <- max(data_sum + data_unc, na.rm = TRUE)
  # Axis settings
  x_axis <- seq(15, 345, by = 30)
  # Plot daily sums
  par(mai = c(0.7, 0.7, 0.7, 0.4)) # set margin
  if (!sum(!is.na(data_sum)) == 0 && missing_unc == 0) {
    # Plot
    plot(
      data_sum ~ data_doy,
      type = "n", ylim = c(y_min, y_max),
      axes = FALSE, xlab = "", ylab = "", main = year
    )
    mtext(unit, 2, 2.2)
    if (var_unc != "none") {
      tb <- !is.na(data_unc) # polygons sensitive to NAs
      polygon(
        c(data_doy[tb], rev(data_doy[tb])),
        c(data_sum[tb] + data_unc[tb], rev(data_sum[tb] - data_unc[tb])),
        col = "dark grey", border = NA
      )
    }
    abline(h = 0, col = "grey")
    lines(data_sum, lty = "solid", lwd = 1, col = "dark green")
    points(data_sum, pch = 20, cex = 0.7, col = "dark green")
    axis(
      1,
      at = x_axis, cex.axis = 1.0, labels = month.abb,
      col.axis = "dark violet"
    )
    axis(2, cex.axis = 1.0)
    box()
  } else {
    # Plot empty box
    plot(
      rep(0, length(data_sum)) ~ data_doy,
      type = "n", axes = FALSE, xlab = "",
      ylab = "", main = year
    )
    axis(
      1,
      at = x_axis, cex.axis = 1.0, labels = month.abb,
      col.axis = "dark violet"
    )
    box()
    if (missing_unc != 0) {
      warning(
        "Uncertainty estimates missing for ", missing_unc,
        " data points of ", var,
        " in year: ", year, ". This will cause an empty plot."
      )
    } else {
      warning("Missing data in year: ", year, ".")
    }
  }
}

plot_windrose <- function(data, ws = ws, wd = wd, ws_res = 0.5, wd_res = 30,
                          ws_min = 0, ws_max = 4, ws_seq = NULL,
                          palette = "Spectral", countmax = NA) {
  ws <- data %>% pull(!!rlang::enquo(ws))
  wd <- data %>% pull(!!rlang::enquo(wd))

  # wd_bin <- ggplot2::cut_number(wd, 360 / wd_res)
  # ws_bin <- ggplot2::cut_interval(ws, length = ws_res)

  # Tidy up input data ----
  n.in <- NROW(data)
  dnu <- (is.na(ws) | is.na(wd))
  ws[dnu] <- NA
  wd[dnu] <- NA

  # figure out the wind speed bins
  if (missing(ws_seq)) {
    ws_seq <- seq(ws_min, ws_max, ws_res)
  } else cat("Using custom speed bins. \n")

  # get some information about the number of bins, etc.
  n_ws_seq <- length(ws_seq)
  n_colors <- n_ws_seq - 1

  # create the color map
  colors <- grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(
      min(max(3, n_colors), min(9, n_colors)), palette
    )))(n_colors)

  if (max(ws, na.rm = TRUE) > ws_max) {
    ws_breaks <- c(ws_seq, max(ws, na.rm = TRUE))
    ws_labs <- c(
      paste(c(ws_seq[1:n_ws_seq - 1]), "-", c(ws_seq[2:n_ws_seq])),
      paste(ws_max, "-", round(max(ws, na.rm = TRUE), 1))
    )
    colors <- c(colors, "grey50")
  } else {
    ws_breaks <- ws_seq
    ws_labs <- paste(c(ws_seq[1:n_ws_seq - 1]), "-", c(ws_seq[2:n_ws_seq]))
  }

  ws_bin <- cut(
    x = ws,
    breaks = ws_breaks,
    labels = ws_labs,
    ordered_result = TRUE
  )

  # figure out the wind direction bins
  wd_breaks <- c(
    -wd_res / 2, seq(wd_res / 2, 360 - wd_res / 2, by = wd_res),
    360 + wd_res / 2
  )
  wd_labs <- c(
    paste(360 - wd_res / 2, "-", wd_res / 2),
    paste(
      seq(wd_res / 2, 360 - 3 * wd_res / 2, by = wd_res), "-",
      seq(3 * wd_res / 2, 360 - wd_res / 2, by = wd_res)
    ),
    paste(360 - wd_res / 2, "-", wd_res / 2)
  )
  # assign each wind direction to a bin
  wd_bin <- cut(wd, breaks = wd_breaks, ordered_result = TRUE)
  levels(wd_bin) <- wd_labs

  # deal with change in ordering introduced somewhere around version 2.2
  if (packageVersion("ggplot2") > "2.2") {
    ws_bin <- factor(ws_bin, levels = rev(levels(ws_bin)))
    colors <- rev(colors)
  }

  # Put everything in a data frame
  df <- data.frame(wd_bin = wd_bin, ws_bin = ws_bin)

  # create the plot
  out <- df %>%
    tidyr::drop_na() %>%
    ggplot2::ggplot(ggplot2::aes(wd_bin, fill = ws_bin)) +
    ggplot2::geom_bar() +
    ggplot2::scale_x_discrete(drop = FALSE, labels = waiver()) +
    ggplot2::coord_polar(start = -((wd_res / 2) / 360) * 2 * pi) +
    ggplot2::scale_fill_brewer(
      name = "Wind speed (m/s)", palette = "Spectral", drop = FALSE
    ) +
    ggplot2::guides(fill = ggplot2::guide_legend(reverse = TRUE)) +
    ggplot2::labs(x = NULL, y = NULL) +
    ggplot2::theme_minimal() +
    ggplot2::theme(
      legend.position = "top",
      axis.text.y = ggplot2::element_blank(),
      panel.grid.major = ggplot2::element_line(colour = "grey80")
    )

  # adjust axes if required
  if (!is.na(countmax)) {
    out <- out + ggplot2::ylim(c(0, countmax))
  }

  out
}

plot_footprint <- function(x, raster = TRUE, contour = TRUE,
                           mask_zero = FALSE) {
  # Hacky solution to list depth issue when fp is in a data frame
  #depth <- purrr::vec_depth(x)
  if (fp_is_empty(x)) stop("No available footprint data.")
  class <- class(x)
  # TODO: give this methods for raster and polygon objects
  if (class == "RasterLayer") {
    fp <- raster::rasterToPoints(x) %>% data.frame()
    fp <- fp %>% dplyr::mutate(z = layer)
    if (mask_zero) fp$z <- replace(fp$z, dplyr::near(fp$z, 0), NA)
  } else {
    percents <- percent_footprint(x)
    matrix <- fp_tidy(x)
    fp <- dplyr::full_join(percents, matrix)
  }

  fp %>%
    ggplot2::ggplot(ggplot2::aes(x, y, z = z)) +
    {if (raster) {
      ggplot2::geom_raster(ggplot2::aes(fill = z), interpolate = TRUE)
    }} +
    {if (contour) {
      ggplot2::stat_contour(
        geom = "contour", color = "black", breaks = c(30, 50, 70, 80, 90, 95)
      )
    }} +
    ggplot2::scale_color_distiller(palette = "RdYlBu", direction = -1) +
    ggplot2::scale_fill_distiller(palette = "RdYlBu", direction = -1) +
    ggplot2::scale_y_continuous(expand = c(0, 0)) +
    ggplot2::scale_x_continuous(expand = c(0, 0)) +
    ggplot2::coord_fixed() +
    #ggplot2::coord_cartesian(xlim = c(-75, 75), ylim = c(-75, 75)) +
    ggplot2::theme_bw()
}

# ==============================================================================
#' Title
#'
#' @param data
#' @param fluxes
#' @param path
#' @param tag
#'
#' @return
#' @export
#'
#' @examples
report_flux <- function(data, fluxes, path = paths$prelim_reports, tag_out) {
  pdf(
    paste0(path, "prelim_report_fluxes_", tag_out, ".pdf"),
    paper = "letter", width = 8.25, height = 10.75, pointsize = 11
  )
  for (i in seq_along(fluxes)) {
    cat(paste0(fluxes[i]), "... ")
    gridExtra::grid.arrange(
      plot_timeseries(
        data, !!as.name(fluxes[i]), means = "daily", n = 7, sample = TRUE
      ) +
        ggplot2::ggtitle(
          paste0(get_full_names()[[fluxes[i]]], " (", fluxes[i], ")"),
          subtitle = "Full time series"
        ) +
        ggplot2::theme(plot.title = ggplot2::element_text(
          face = "bold", hjust = 0.5, size = 14
        )),
      plot_heatmap(data, !!as.name(fluxes[i])) +
        ggplot2::labs(subtitle = "Fingerprint heatmap"),
      plot_daily_means(
        data, !!as.name(fluxes[i]), n = 7, title = TRUE
      ),
      ncol = 1, heights = c(0.3, 0.3, 0.4)
    )
    gridExtra::grid.arrange(
      plot_diurnal(data, !!as.name(fluxes[i]), groups = "month", light = ppfd) +
        ggplot2::ggtitle(
          paste0(get_full_names()[[fluxes[i]]], " (", fluxes[i], ")"),
          subtitle = "Monthly diurnal cycles"
        ) +
        ggplot2::theme(plot.title = ggplot2::element_text(
          face = "bold", hjust = 0.5, size = 14
        )),
      plot_wd_stats(
        data, !!as.name(fluxes[i]), groups = "daynight", sample = TRUE
      ) +
        ggplot2::labs(subtitle = "Wind direction dependency"),
      plot_covariates(
        data, !!as.name(fluxes[i]), filter_qc = TRUE, sample = TRUE
      ) +
        ggplot2::labs(subtitle = "Strongest environmental predictors"),
      ncol = 1, heights = c(0.4, 0.3, 0.3)
    )
  }
  quiet(dev.off())
  cat("done")
}

# ==============================================================================
#' Title
#'
#' @param data
#' @param vars
#' @param path
#' @param tag
#'
#' @return
#' @export
#'
#' @examples
report_vars <- function(data, vars, path = paths$prelim_reports, tag_out,
                        means = "daily", remove_outliers = TRUE, z = 7,
                        sample = TRUE, show_light = TRUE, light = NULL) {
  pdf(
    paste0(path, tag_out, ".pdf"),
    paper = "letter", width = 8.25, height = 10.75, pointsize = 11
  )
  for (i in seq_along(vars)) {
    cat(paste0(vars[i]), "... ")
    gridExtra::grid.arrange(
      plot_timeseries(
        data, !!as.name(vars[i]), means = means,
        remove_outliers = remove_outliers, z = z, sample = sample
      ) +
        ggplot2::ggtitle(
          paste0(get_full_names()[[vars[i]]], " (", vars[i], ")"),
          subtitle = "Full time series"
        ) +
        ggplot2::theme(plot.title = ggplot2::element_text(
          face = "bold", hjust = 0.5, size = 14
        )),
      gridExtra::arrangeGrob(
        plot_diurnal(
          data, !!as.name(vars[i]), remove_outliers = remove_outliers, z = z,
          show_light = show_light, light = light
        ) +
          ggplot2::labs(subtitle = "Average diurnal cycle"),
        plot_distribution(
          data, !!as.name(vars[i]), groups = "daynight",
          remove_outliers = remove_outliers, z = z
        ) +
          ggplot2::labs(subtitle = "Day/night distributions"),
        ncol = 2
      ),
      plot_wd_stats(
        data, !!as.name(vars[i]), remove_outliers = remove_outliers, z = z,
        sample = sample
      ) +
        ggplot2::labs(subtitle = "Wind direction dependency"),
      ncol = 1, heights = c(0.4, 0.3, 0.3)
    )
  }
  quiet(dev.off())
  cat("done")
}

report_site <- function(data, path, tag_out) {
  pdf(
    paste0(path, "prelim_report_site_", tag_out, ".pdf"),
    paper = "letter", width = 8.25, height = 10.75, pointsize = 11
  )
  for (i in seq_along(vars)) {
    cat(paste0(vars[i]), "... ")
    gridExtra::grid.arrange(
      "footprint",
      "wind rose",
      plot3,
      ncol = 1, heights = c(0.4, 0.3, 0.3)
    )
  }
  quiet(dev.off())
  cat("done")
}