Dev

R/GetVarFromSpectra.R

@@ -41,8 +41,9 @@ GetVarFromSpectra <- function(spec,f,dfreq=NULL,df.log = 0,bw = 3)
 
     if (df.log > 0) spec <- LogSmooth(spec,removeLast = 0,df.log = df.log)
 
-     vars <- mean(SpecInterpolate(spec,newFreq)$spec)*diff(f)*2
-     dof <- mean(SpecInterpolate(spec,newFreq)$dof)
+     spec.int <- SpecInterpolate(spec, newFreq)
+     vars <- mean(spec.int$spec)*diff(f)*2
+     dof <- mean(spec.int$dof)
 
      ## Estimate the DOF by calculating how many independent
      ## spectral estimates contribute to the calculated mean value

R/gg_spec.R

@@ -6,15 +6,17 @@
 #' @param conf Plot shaded confidence interval if it exists in the spec object
 #' @param alpha.line Alpha level for the spectra line(s)
 #' @param alpha.ribbon Alpha level for the confidence region(s)
-#' @param colour Colour or name of variable to which to map to colour
-#' @param min.colours minimum number of spectra before starting to colour them
+#' @param colour Name of variable map to colour, unquoted
+#' @param linetype Name of variable map to linetype, unquoted
+#' @param group Name of variable to map to group, unquoted
+#' @param min.colours Minimum number of spectra before starting to colour them
 #' separately
-#' @param removeFirst,removeLast remove first or last "n" values on the low or
-#'   high frequency side respectively. Will be unpredictable when used with a
-#'   list of spectra, or spec_df object with multiple spectra with different
-#'   frequency axes.
-#' @param force.CI Force the plotting of confidence regions when the total number
-#' of frequencies exceeds 10000. Defaults to FALSE
+#' @param removeFirst,removeLast Remove first or last "n" values on the low or
+#'   high frequency side respectively. Operates on a per group basis.
+#' @param force.lims,force.CI Force the plotting of confidence regions when the total number
+#' of frequencies exceeds 10000. force.CI is deprecated, use force.lims. Defaults to FALSE
+#' @param quantiles Plot uncertainty quantiles if they exist in the spec object
+#' @param time_unit Optional string giving the time unit for axes labels, e.g. "years"
 #'
 #' @return a ggplot object
 #' @family functions to plot power spectra
@@ -24,13 +26,13 @@
 #' library(PaleoSpec)
 #' N <- 1e03
 #' beta <- 1
-#' alpha = 0.1
+#' alpha <- 0.1
 #'
 #' ts1 <- SimPLS(N = N, b = beta, a = alpha)
 #' ts2 <- SimPLS(N = N, b = beta, a = alpha)
-#' sp1 <- SpecMTM(ts1)
+#' sp1 <- SpecMTM(ts1, deltat = 1)
 #' sp1 <- AddConfInterval(sp1)
-#' sp2 <- SpecMTM(ts2)
+#' sp2 <- SpecMTM(ts2, deltat = 1)
 #'
 #' # plot single spectrum
 #' p <- gg_spec(sp1, spec_id = "df1")
@@ -59,47 +61,50 @@
 #' sp2 <- LogSmooth(sp1)
 #' p <- gg_spec(sp2, p)
 #' p
-gg_spec <- function(x, gg = NULL, conf = TRUE,
+gg_spec <- function(x, gg = NULL,
+                    conf = TRUE,
                     spec_id = NULL,
                     colour = spec_id,
+                    group = spec_id,
+                    linetype = NULL,
                     alpha.line = 1,
-                    alpha.ribbon = 0.3,
+                    alpha.ribbon = c(0.166, 0.333),
                     removeFirst = 0, removeLast = 0,
-                    force.CI = FALSE,
-                    min.colours = 2) {
+                    min.colours = 2,
+                    force.lims = FALSE,
+                    force.CI = NULL,
+                    quantiles = FALSE,
+                    time_unit = NULL) {
+  if (!missing("force.CI")) {
+
+    warning("the force.CI argument is deprecated please use force.lims")
+    force.lims <- force.CI
+  }
 
   gg_installed <- requireNamespace("ggplot2", quietly = TRUE)
 
-  if (gg_installed == FALSE){
+  if (gg_installed == FALSE) {
     stop("package ggplot2 is required to use gg_spec(). To install ggplot2, run install.packages(\"ggplot2\") from the console")
   }
 
-  if (class(x)[1] != "spec_df" & "list" %in% class(x) == FALSE){
+  if (class(x)[1] != "spec_df" & "list" %in% class(x) == FALSE) {
     x <- list(x)
     names(x) <- spec_id
   }
 
-  if (class(x)[1] != "spec_df"){
+  if (class(x)[1] != "spec_df") {
     df <- Spec2DF(x)
   } else {
     df <- x
   }
 
 
-  ## Remove first and or last freqs on a per spec_id basis
-
-  if (exists("spec_id", df) == FALSE){
-    df$spec_id <- 1
-  }
-
   if (is.numeric(df$spec_id)){
     df$spec_id <- as.character(df$spec_id)
   }
 
 
   if (removeFirst > 0) {
-    #df <- df[rank(df$freq) > removeFirst, ]
-
     # Convert 'df' to a data.table
     data.table::setDT(df)
 
@@ -111,7 +116,6 @@ gg_spec <- function(x, gg = NULL, conf = TRUE,
   }
 
   if (removeLast > 0) {
-    #df <- df[rank(-df$freq) > removeLast,]
     # Convert 'df' to a data.table
     data.table::setDT(df)
 
@@ -120,68 +124,107 @@ gg_spec <- function(x, gg = NULL, conf = TRUE,
 
     # return to data.frame
     df <- as.data.frame(df)
+
+  }
+
+  if (is.null(time_unit) == FALSE) {
+    lab_timescale <- paste0("Timescale [", time_unit, "]")
+    lab_freq <- paste0("Frequency [1/", time_unit, "]")
+  } else {
+    lab_timescale <- "Timescale"
+    lab_freq <- "Frequency"
   }
 
 
   if (is.null(gg)) {
-    p <- ggplot2::ggplot(data = df)
-  } else {
+    p <- ggplot2::ggplot(data = df, ggplot2::aes(group = {{ group }}))
+    } else {
     p <- gg
   }
 
   # rename to PSD so that y axis label can be overwritten later
   df$PSD <- df$spec
   df$Frequency <- df$freq
 
-  if (conf == TRUE & exists("lim.1", df)){
+  df <- as.data.frame(df)
 
-    if (nrow(df) > 1e04 & force.CI == FALSE){
+  if (conf == TRUE & exists("lim.1", df)) {
+    if (nrow(df) > 1e04 & force.lims == FALSE) {
       warning("geom_ribbon is very slow when the number of points > 1e04, skipping the confidence region")
     } else {
+    p <- p +
+      ggplot2::geom_ribbon(
+          data = df, ggplot2::aes(
+            x = Frequency, ymin = lim.2,
+            ymax = lim.1, fill = {{ colour }}
+          ),
+          alpha = alpha.ribbon[1], colour = NA
+        )
+    }
+    }
 
+
+  if (quantiles == TRUE & exists("X2.5.", df)) {
+    if (nrow(df) > 1e04 & force.lims == FALSE) {
+      warning("geom_ribbon is very slow when the number of points > 1e04, skipping the confidence region")
+    } else {
       p <- p +
-        ggplot2::geom_ribbon(data = df, ggplot2::aes(x = Frequency, ymin = lim.2,
-                                                     group = spec_id,
-                                                     ymax = lim.1, fill = {{ colour }}),
-                             alpha = alpha.ribbon, colour = NA)
+        ggplot2::geom_ribbon(
+          data = df, ggplot2::aes(
+            x = Frequency, ymin = `X2.5.`,
+            ymax = `X97.5.`, fill = {{ colour }}),
+          alpha = alpha.ribbon[1], colour = NA
+        ) +
+        ggplot2::geom_ribbon(
+          data = df, ggplot2::aes(
+            x = Frequency, ymin = `X15.9.`,
+            ymax = `X84.1.`, fill = {{ colour }}
+          ),
+          alpha = alpha.ribbon[2], colour = NA
+        )
     }
   }
 
-  p <- p + ggplot2::geom_line(data = df, ggplot2::aes(x = Frequency, y = PSD,
-                                                      group = spec_id,
-                                                      colour = {{ colour }}),
-                              alpha = alpha.line
+  p <- p + ggplot2::geom_line(
+    data = df, ggplot2::aes(
+      x = Frequency, y = PSD,
+      linetype = {{ linetype }},
+      colour = {{ colour }}
+    ),
+    alpha = alpha.line
   ) +
-    ggplot2::scale_x_continuous(trans = "log10",
-                                sec.axis = ggplot2::sec_axis(~ 1/., "Timescale")) +
+    ggplot2::scale_x_continuous(
+      name = lab_freq, trans = "log10",
+      sec.axis = ggplot2::sec_axis(~ 1 / ., lab_timescale, labels = scales::comma)
+    ) +
     ggplot2::scale_y_continuous(trans = "log10") +
     ggplot2::annotation_logticks(sides = "tlb") +
     ggplot2::theme_bw() +
     ggplot2::scale_colour_brewer("",
-                                 type = "qual",
-                                 palette = "Dark2",
-                                 aesthetics = c("colour", "fill")) +
+      type = "qual",
+      palette = "Dark2",
+      aesthetics = c("colour", "fill")
+    ) +
     ggplot2::theme(panel.grid.minor = ggplot2::element_blank()) +
-    ggplot2::scale_alpha()
+    ggplot2::scale_alpha() +
+    ggplot2::labs(linetype = "")
 
   g <- ggplot2::ggplot_build(p)
   colrs <- unlist(unique(sapply(g$data, function(x) unique(x["colour"])$colour)))
   colrs <- colrs[is.na(colrs) == FALSE]
   ncolrs <- length(colrs)
 
-  if (ncolrs <= min.colours){
 
-    p <- p + ggplot2::scale_colour_manual("", values = "black",
-                                          aesthetics = c("colour", "fill"))
+  if (ncolrs <= min.colours) {
+    p <- p + ggplot2::scale_colour_manual("",
+      values = "black",
+      aesthetics = c("colour", "fill")
+    )
 
-    if (is.null({{ colour }})){
+    if (is.null({{ colour }}) & is.null({{ linetype }})) {
       p <- p + ggplot2::theme(legend.position = "none")
     }
-
   }
 
   p
 }
-
-
-