R/specta_data_ploting.R

In resourcecode: Access to the 'RESOURCECODE' Hindcast Database

#' Plot a wave density 2D spectrum
#'
#' @param spec the spectral data, as an output from `get_2Dspectrum`
#' @param time the time to plot. Either an integer or the date.
#' @param normalize Should the spectrum be normalized to have maimum 1 before ploting
#' @param trim removes the values of the spectral density lower than this value
#' @param cut_off cut-off frequency above which the spectrum is not plotted
#' @param ... currently unused
#'
#' @return a ggplot object
#' @export
#'
#' @examplesIf requireNamespace("resourcecodedata", quietly = TRUE)
#' spec <- get_2d_spectrum("SEMREVO", start = "1994-01-01", end = "1994-01-31")
#' plot_2d_specta(spec, 1)
#' @importFrom ggplot2 ggplot geom_rect scale_x_continuous scale_y_continuous
#'                     theme_bw coord_polar scale_color_distiller scale_fill_distiller
#'                     labs expansion
plot_2d_specta <- function(
  spec,
  time = 1L,
  normalize = TRUE,
  trim = 0.01,
  cut_off = 0.4,
  ...
) {
  if (is.character(time)) {
    time <- as.POSIXct(time, tz = "UTC")
  }

  if (inherits(time, "POSIXct")) {
    time <- which(time == spec$forcings$time)
  }

  df_freq <- tidyr::expand_grid(frequency1 = 1:36, dir = spec$dir)

  df_freq$frequency2 <- spec$frequency2[df_freq$frequency1]
  df_freq$frequency1 <- spec$frequency1[df_freq$frequency1]

  df <- tibble::tibble(df_freq, ef = c(spec$efth[, , time]))

  legend_text <- "Power spectrum\n(m^2.s)"

  if (normalize) {
    ef_max <- max(df$ef)
    df$ef <- df$ef / ef_max

    legend_text <- "Normalized\nPower spectrum"
  }

  if (!is.null(trim)) {
    # Trim
    df$ef[df$ef <= trim] <- NA
  }

  df$dir <- (df$dir + 180) %% 360

  ggplot(
    df,
    aes(
      xmin = .data$frequency1,
      xmax = .data$frequency2,
      ymin = dir - 5,
      ymax = dir + 5,
      fill = .data$ef,
      col = .data$ef
    )
  ) +
    geom_rect() +
    scale_x_continuous(
      name = "Frequency (Hz)",
      expand = expansion(),
      limits = c(0, cut_off)
    ) +
    scale_y_continuous(
      name = "Direction from (\u00b0)",
      expand = expansion(),
      breaks = c(0, 90, 180, 270),
      labels = c("N", "E", "S", "W"),
      minor_breaks = seq(from = 0, to = 360, by = 30)
    ) +
    coord_polar(theta = "y", start = -5 * pi / 180) +
    scale_color_distiller(
      palette = "PuBu",
      direction = 1,
      name = legend_text,
      na.value = "transparent"
    ) +
    scale_fill_distiller(
      palette = "PuBu",
      direction = 1,
      name = legend_text,
      na.value = "transparent"
    ) +
    labs(
      title = paste(
        "Directional Wave energy spectrum at location",
        spec$station
      ),
      subtitle = format(spec$forcings$time[time], format = "%Y-%m-%d %H:%M"),
      caption = "Source: Resourcecode hindcast database\nresourcecode.ifremer.fr"
    ) +
    theme_bw() +
    theme(legend.position = "bottom")
}