R/visualization.R

In grahamstewart12/fluxtools:

## 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_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")
}