R/plot.biwavelet.R

In biwavelet: Conduct Univariate and Bivariate Wavelet Analyses

#' Plot \code{biwavelet} objects
#'
#' Plot \code{biwavelet} objects such as the cwt, cross-wavelet and wavelet
#' coherence.
#'
#' @param x \code{biwavelet} object generated by \code{\link{wt}},
#'   \code{\link{xwt}}, or \code{\link{wtc}}.
#' @param ncol Number of colors to use.
#' @param fill.cols Vector of fill colors to be used. Users can specify color
#'   vectors using \code{\link{colorRampPalette}} or
#'   \code{\link[RColorBrewer]{brewer.pal}} from package
#'   \code{\link[RColorBrewer]{RColorBrewer}}. Value \code{NULL} generates
#'   MATLAB's jet color palette.
#' @param xlab X-label of the figure.
#' @param ylab Y-label of the figure.
#' @param tol Tolerance level for significance contours. Sigificance contours
#'   will be drawn around all regions of the spectrum where
#'   \code{spectrum/percentile >= tol}. If strict \eqn{i^{th}} percentile
#'   regions are desired, then \code{tol} must be set to 1.
#' @param plot.cb Plot color bar if \code{TRUE}.
#' @param plot.phase Plot phases with black arrows.
#' @param type Type of plot to create. Can be \code{power} to plot the power,
#'   \code{power.corr} to plot the bias-corrected power, \code{power.norm} to
#'   plot the power normalized by the variance, \code{power.corr.norm} to plot
#'   the bias-corrected power normalized by the variance, \code{wavelet} to plot
#'   the wavelet coefficients, or \code{phase} to plot the phase.
#' @param plot.coi Plot cone of influence (COI) as a semi-transparent polygon if
#'   \code{TRUE}. Areas that fall within the polygon can be affected by edge
#'   effects.
#' @param lwd.coi Line width of COI.
#' @param col.coi Color of COI.
#' @param lty.coi Line type of COI. Value 1 is for solide lines.
#' @param alpha.coi Transparency of COI. Range is 0 (full transparency) to 1 (no
#'   transparency).
#' @param plot.sig Plot contours for significance if \code{TRUE}.
#' @param lwd.sig Line width of significance contours.
#' @param col.sig Color of significance contours.
#' @param lty.sig Line type of significance contours.
#' @param bw plot in black and white if \code{TRUE}.
#' @param legend.loc Legend location coordinates as defined by
#'   \code{\link[fields]{image.plot}}.
#' @param legend.horiz Plot a horizontal legend if \code{TRUE}.
#' @param arrow.len Size of the arrows. Default is based on plotting region.
#' @param arrow.lwd Width/thickness of arrows.
#' @param arrow.cutoff Cutoff value for plotting phase arrows. Phase arrows will
#'   be be plotted in regions where the significance of the zvalues exceeds
#'   \code{arrow.cutoff} for \code{\link{wt}} and \code{\link{xwt}} objects. For
#'   \code{\link{pwtc}} and \code{\link{wtc}} objects, phase arrows will be
#'   plotted in regions where the \code{rsq} value exceeds \code{arrow.cutoff}.
#'   If the object being plotted does not have a significance field, regions
#'   whose z-values exceed the \code{arrow.cutoff} quantile will be plotted.
#' @param arrow.col Color of arrows.
#' @param xlim The x limits.
#' @param ylim The y limits.
#' @param zlim The z limits.
#' @param xaxt Add x-axis? Use \code{n} for none.
#' @param yaxt Add y-axis? Use \code{n} for none.
#' @param form Format to use to display dates on the x-axis.
#' See \code{\link{Date}} for other valid formats.
#' @param \dots Other parameters.
#'
#' @details
#' Arrows pointing to the right mean that \code{x} and \code{y} are in phase.
#'
#' Arrows pointing to the left mean that \code{x} and \code{y} are in anti-phase.
#'
#' Arrows pointing up mean that \code{x} leads \code{y} by \eqn{\pi/2}.
#'
#' Arrows pointing down mean that \code{y} leads \code{x} by \eqn{\pi/2}.
#'
#' @references
#' Cazelles, B., M. Chavez, D. Berteaux, F. Menard, J. O. Vik, S. Jenouvrier,
#' and N. C. Stenseth. 2008. Wavelet analysis of ecological time series.
#' \emph{Oecologia} 156:287-304.
#'
#' Grinsted, A., J. C. Moore, and S. Jevrejeva. 2004. Application of the cross
#' wavelet transform and wavelet coherence to geophysical time series.
#' \emph{Nonlinear Processes in Geophysics} 11:561-566.
#'
#' Torrence, C., and G. P. Compo. 1998. A Practical Guide to Wavelet Analysis.
#' \emph{Bulletin of the American Meteorological Society} 79:61-78.
#'
#' Liu, Y., X. San Liang, and R. H. Weisberg. 2007. Rectification of the Bias in
#' the Wavelet Power Spectrum. \emph{Journal of Atmospheric and Oceanic
#' Technology} 24:2093-2102.
#'
#' @author Tarik C. Gouhier (tarik.gouhier@@gmail.com)
#' Code based on WTC MATLAB package written by Aslak Grinsted.
#'
#' @examples
#' t1 <- cbind(1:100, rnorm(100))
#' t2 <- cbind(1:100, rnorm(100))
#'
#' # Continuous wavelet transform
#' wt.t1 <- wt(t1)
#'
#' # Plot power
#' # Make room to the right for the color bar
#' par(oma = c(0, 0, 0, 1), mar = c(5, 4, 4, 5) + 0.1)
#' plot(wt.t1, plot.cb = TRUE, plot.phase = FALSE)
#'
#' # Cross-wavelet transform
#' xwt.t1t2 <- xwt(t1, t2)
#'
#' # Plot cross-wavelet
#' par(oma = c(0, 0, 0, 1), mar = c(5, 4, 4, 5) + 0.1)
#' plot(xwt.t1t2, plot.cb = TRUE)
#'
#' # Example of bias-correction
#' t1 <- sin(seq(0, 2 * 5 * pi, length.out = 1000))
#' t2 <- sin(seq(0, 2 * 15 * pi, length.out = 1000))
#' t3 <- sin(seq(0, 2 * 40 * pi, length.out = 1000))
#'
#' # This aggregate time series should have the same power
#' # at three distinct periods
#' s <- t1 + t2 + t3
#'
#' # Compare plots to see bias-effect on CWT:
#' # biased power spectrum artificially
#' # reduces the power of higher-frequency fluctuations.
#' wt1 <- wt(cbind(1:1000, s))
#' par(mfrow = c(1,2))
#' plot(wt1, type = "power.corr.norm", main = "Bias-corrected")
#' plot(wt1, type = "power.norm", main = "Biased")
#'
#' # Compare plots to see bias-effect on XWT:
#' # biased power spectrum artificially
#' # reduces the power of higher-frequency fluctuations.
#' x1 <- xwt(cbind(1:1000, s), cbind(1:1000, s))
#' par(mfrow = c(1,2))
#'
#' plot(x1, type = "power.corr.norm", main = "Bias-corrected")
#' plot(x1, type = "power.norm", main = "Biased")
#'
#' @importFrom fields image.plot
#' @importFrom grDevices adjustcolor colorRampPalette
#' @importFrom graphics axTicks axis box contour image par polygon
#' @export
plot.biwavelet <- function(x, ncol = 64, fill.cols = NULL,
                           xlab = "Time", ylab = "Period",
                           tol = 1, plot.cb = FALSE, plot.phase = FALSE,
                           type = "power.corr.norm",
                           plot.coi = TRUE, lwd.coi = 1, col.coi = "white",
                           lty.coi = 1, alpha.coi = 0.5, plot.sig = TRUE,
                           lwd.sig = 4, col.sig = "black", lty.sig = 1,
                           bw = FALSE,
                           legend.loc = NULL,
                           legend.horiz = FALSE,
                           arrow.len = min(par()$pin[2] / 30,
                                           par()$pin[1] / 40),
                           arrow.lwd = arrow.len * 0.3,
                           arrow.cutoff = 0.8,
                           arrow.col = "black",
                           xlim = NULL, ylim = NULL, zlim = NULL,
                           xaxt = "s", yaxt = "s", form = "%Y", ...) {

  if (is.null(fill.cols)) {
    if (bw) {
      fill.cols <- c("black", "white")
    } else {
      fill.cols <- c("#00007F", "blue", "#007FFF",
                     "cyan","#7FFF7F", "yellow",
                     "#FF7F00", "red", "#7F0000")
    }
  }

  col.pal <- colorRampPalette(fill.cols)
  fill.colors <- col.pal(ncol)

  types <- c("power.corr.norm", "power.corr", "power.norm",
             "power", "wavelet", "phase", "timing.err")

  type <- match.arg(tolower(type), types)

  if (type == "power.corr" | type == "power.corr.norm") {
    if (x$type == "wtc" | x$type == "xwt") {
      x$power <- x$power.corr
      x$wave <- x$wave.corr
    } else {
      x$power <- x$power.corr
    }
  }

  if (type == "power.norm" | type == "power.corr.norm") {
    if (x$type == "xwt") {
      zvals <- log2(x$power) / (x$d1.sigma * x$d2.sigma)

      if (is.null(zlim)) {
        zlim <- range(c(-1, 1) * max(zvals))
      }

      zvals[zvals < zlim[1]] <- zlim[1]
      locs <- pretty(range(zlim), n = 5)
      leg.lab <- 2 ^ locs

    } else if (x$type == "wtc" | x$type == "pwtc") {
      zvals <- x$rsq
      zvals[!is.finite(zvals)] <- NA
      if (is.null(zlim)) {
        zlim <- range(zvals, na.rm = TRUE)
      }
      zvals[zvals < zlim[1]] <- zlim[1]
      locs <- pretty(range(zlim), n = 5)
      leg.lab <- locs
    } else {
      zvals <- log2(abs(x$power / x$sigma2))
      if (is.null(zlim)) {
        zlim <- range(c(-1, 1) * max(zvals))
      }
      zvals[zvals < zlim[1]] <- zlim[1]
      locs <- pretty(range(zlim), n = 5)
      leg.lab <- 2 ^ locs
    }
  } else if (type == "power" | type == "power.corr") {
    zvals <- log2(x$power)
    if (is.null(zlim)) {
      zlim <- range( c(-1, 1) * max(zvals) )
    }
    zvals[zvals < zlim[1]] <- zlim[1]
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- 2 ^ locs
  } else if (type == "wavelet") {
    zvals <- (Re(x$wave))
    if (is.null(zlim)) {
      zlim <- range(zvals)
    }
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- locs
  } else if (type == "phase") {
    zvals <- x$phase
    if (is.null(zlim)) {
      zlim <- c(-pi, pi)
    }
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- locs
  } else if (type == "timing.err") {
    zvals <- x$timing.err
    if (is.null(zlim)) {
      zlim <- range(zvals)
    }
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- locs
  }


  if (is.null(xlim)) {
    xlim <- range(x$t)
  }

  yvals <- log2(x$period)

  if (is.null(ylim)) {
    ylim <- range(yvals)
  } else {
    ylim <- log2(ylim)
  }

  image(x$t,
        yvals,
        t(zvals),
        zlim = zlim,
        xlim = xlim,
        ylim = rev(ylim),
        xlab = xlab,
        ylab = ylab,
        yaxt = "n",
        xaxt = "n",
        col = fill.colors, ...)

  box()
  if (class(x$xaxis)[1] == "Date" | class(x$xaxis)[1] == "POSIXct") {
    if (xaxt != "n") {
      xlocs <- pretty(x$t) + 1
      axis(side = 1, at = xlocs, labels = format(x$xaxis[xlocs], form, ...))
    }
  } else {
    if (xaxt != "n") {
      xlocs <- axTicks(1)
      axis(side = 1, at = xlocs, ...)
    }
  }

  if (yaxt != "n") {
    axis.locs <- axTicks(2)
    yticklab <- format(2 ^ axis.locs) #, dig = 1)
    axis(2, at = axis.locs, labels = yticklab, ...)
  }

  # COI
  if (plot.coi) {
    # polygon(x = c(x$t, rev(x$t)), lty = lty.coi, lwd = lwd.coi,
    #         y = c(log2(x$coi),
    #               rep(max(log2(x$coi), na.rm = TRUE), length(x$coi))),
    #         col = adjustcolor(col.coi, alpha.f = alpha.coi), border = col.coi)
    polygon(x = c(x$t, rev(x$t)), lty = lty.coi, lwd = lwd.coi,
            y = c(log2(x$coi), rep(max(c(log2(x$coi), x$period), na.rm = TRUE),
                                   length(x$coi))),
            col = adjustcolor(col.coi, alpha.f = alpha.coi), border = col.coi)
  }

  # sig.level contour
  if (plot.sig & length(x$signif) > 1) {
    if (x$type %in% c("wt", "xwt")) {
      contour(x$t, yvals, t(x$signif), level = tol, col = col.sig,
              lwd = lwd.sig, add = TRUE, drawlabels = FALSE)
    } else {
      tmp <- x$rsq / x$signif
      contour(x$t, yvals, t(tmp), level = tol, col = col.sig, lwd = lwd.sig,
              add = TRUE, drawlabels = FALSE)
    }
  }

  # Plot phases
  if (plot.phase) {
    a <- x$phase

    # Remove phases where power is weak
    if (!is.null(x$type)) {
      if (x$type %in% c("wt", "xwt")) {
        locs.phases <- which(x$signif <= arrow.cutoff)
      } else if (x$type %in% c("wtc", "pwtc")) {
        # v <- x$rsq / x$signif
        v <- x$rsq
        locs.phases <- which(v <= arrow.cutoff)
      }
    } else {
      locs.phases <- which(zvals < quantile(zvals, arrow.cutoff))
    }
    a[locs.phases] <- NA

    phase.plot(x$t, log2(x$period), phases = a,
               arrow.len = arrow.len,
               arrow.lwd = arrow.lwd,
               arrow.col = arrow.col)
  }
  box()

  ## Add color bar: this must happen after everything, otherwise chaos ensues!
  if (plot.cb) {
    image.plot(x$t,
               yvals,
               t(zvals),
               zlim = zlim,
               ylim = rev(range(yvals)),
               xlab = xlab,
               ylab = ylab,
               col = fill.colors,
               smallplot = legend.loc,
               horizontal = legend.horiz,
               legend.only = TRUE,
               axis.args =
                 list(at = locs, labels = format(leg.lab, dig = 2), ...),
               xpd = NA)
  }
}