#' 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.
#'
#' @seealso \code{\link{image.plot}}
#'
#' @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)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.