wt.image <- function (WT, my.series = 1, exponent = 1, plot.coi = TRUE, plot.contour = TRUE,
siglvl = 0.1, col.contour = "white", plot.ridge = TRUE, lvl = 0,
col.ridge = "black", color.key = "quantile", n.levels = 100,
color.palette = "rainbow(n.levels, start = 0, end = .7)",
maximum.level = NULL, useRaster = TRUE, max.contour.segments = 250000,
plot.legend = TRUE, legend.params = list(width = 1.2, shrink = 0.9,
mar = 5.1, n.ticks = 6, label.digits = 1, label.format = "f",
lab = NULL, lab.line = 2.5), label.time.axis = TRUE,
show.date = FALSE, date.format = NULL, date.tz = NULL, timelab = NULL,
timetck = 0.02, timetcl = 0.5, spec.time.axis = list(at = NULL,
labels = TRUE, las = 1, hadj = NA, padj = NA), label.period.axis = TRUE,
periodlab = NULL, periodtck = 0.02, periodtcl = 0.5, spec.period.axis = list(at = NULL,
labels = TRUE, las = 1, hadj = NA, padj = NA), main = NULL,
lwd = 2, lwd.axis = 1, graphics.reset = TRUE, verbose = FALSE)
{
if (verbose == T) {
out <- function(...) {
cat(...)
}
}
else {
out <- function(...) {
}
}
default.options = options()
options(max.contour.segments = as.integer(max.contour.segments))
axis.1 <- WT$axis.1
axis.2 <- WT$axis.2
series.data = WT$series
if (exponent <= 0) {
stop("Please use a positive exponent, or return to default setting (exponent=1)!")
}
if (class(WT) == "analyze.wavelet") {
out("Your input object class is 'analyze.wavelet'...\n")
my.series = ifelse(names(series.data)[1] == "date", names(series.data)[2],
names(series.data)[1])
Power = WT$Power^exponent
Power.pval = WT$Power.pval
Ridge = WT$Ridge
}
if (class(WT) == "analyze.coherency") {
out("Your input object class is 'analyze.coherency'...\n")
if (is.numeric(my.series)) {
if (!is.element(my.series, c(1, 2))) {
stop("Please choose either series number 1 or 2!")
}
my.series = ifelse(names(series.data)[1] == "date",
names(series.data)[my.series + 1], names(series.data)[my.series])
}
ind = which(names(series.data) == my.series)
which.series.num = ifelse(names(series.data)[1] == "date",
ind - 1, ind)
if (!is.element(which.series.num, c(1, 2))) {
stop("Your series name is not available, please check!")
}
if (which.series.num == 1) {
Power = WT$Power.x^exponent
Power.pval = WT$Power.x.pval
Ridge = WT$Ridge.x
}
if (which.series.num == 2) {
Power = WT$Power.y^exponent
Power.pval = WT$Power.y.pval
Ridge = WT$Ridge.y
}
}
out(paste("A wavelet power image of your time series '",
my.series, "' will be plotted...", sep = ""), "\n")
if ((!is.null(maximum.level)) & (is.element(color.key, c("quantile",
"q")))) {
warning("\nPlease set color.key = 'i' to make your maximum level specification effective.",
immediate. = TRUE)
}
if (is.element(color.key, c("interval", "i"))) {
if (is.null(maximum.level)) {
maximum.level = max(Power)
}
if (maximum.level < max(Power)) {
stop(paste("... plot can't be produced! Your choice of maximum plot level is smaller than the maximum level observed! Please choose maximum.level larger than ",
max(Power), " or return to default setting (maximum.level = NULL)!",
sep = ""))
}
}
if (is.element(color.key, c("interval", "i"))) {
wavelet.levels = seq(from = 0, to = maximum.level, length.out = n.levels +
1)
}
if (is.element(color.key, c("quantile", "q"))) {
wavelet.levels = quantile(Power, probs = seq(from = 0,
to = 1, length.out = n.levels + 1))
}
key.cols = rev(eval(parse(text = color.palette)))
if (!is.list(legend.params))
legend.params = list()
if (is.null(legend.params$width))
legend.params$width = 1.2
if (is.null(legend.params$shrink))
legend.params$shrink = 0.9
if (is.null(legend.params$mar))
legend.params$mar = ifelse(is.null(legend.params$lab),
5.1, 6.1)
if (is.null(legend.params$n.ticks))
legend.params$n.ticks = 6
if (is.null(legend.params$label.digits))
legend.params$label.digits = 1
if (is.null(legend.params$label.format))
legend.params$label.format = "f"
if (is.null(legend.params$lab.line))
legend.params$lab.line = 2.5
if (is.null(date.format)) {
date.format = WT$date.format
}
if (is.null(date.tz)) {
date.tz = ifelse(is.null(WT$date.tz), "", WT$date.tz)
}
if (!is.list(spec.time.axis))
spec.time.axis = list()
if (is.null(spec.time.axis$at))
spec.time.axis$at = NULL
if (is.null(spec.time.axis$labels))
spec.time.axis$labels = T
if (is.null(spec.time.axis$las))
spec.time.axis$las = 1
if (is.null(spec.time.axis$hadj))
spec.time.axis$hadj = NA
if (is.null(spec.time.axis$padj))
spec.time.axis$padj = NA
time.axis.warning = F
time.axis.warning.na = F
chronology.warning = F
if ((!is.null(spec.time.axis$at)) & (label.time.axis == F))
warning("\nPlease set label.time.axis = TRUE to make time axis specification effective.",
immediate. = TRUE)
if (!is.list(spec.period.axis))
spec.period.axis = list()
if (is.null(spec.period.axis$at))
spec.period.axis$at = NULL
if (is.null(spec.period.axis$labels))
spec.period.axis$labels = T
if (is.null(spec.period.axis$las))
spec.period.axis$las = 1
if (is.null(spec.period.axis$hadj))
spec.period.axis$hadj = NA
if (is.null(spec.period.axis$padj))
spec.period.axis$padj = NA
period.axis.warning = F
period.axis.warning.na = F
if ((!is.null(spec.period.axis$at)) & (label.period.axis ==
F))
warning("\nPlease set label.period.axis = TRUE to make period axis specification effective.",
immediate. = TRUE)
op = par(no.readonly = TRUE)
image.plt = par()$plt
legend.plt = NULL
if (plot.legend == T) {
legend.plt = par()$plt
char.size = par()$cin[1]/par()$din[1]
hoffset = char.size * par()$mar[4]
legend.width = char.size * legend.params$width
legend.mar = char.size * legend.params$mar
legend.plt[2] = 1 - legend.mar
legend.plt[1] = legend.plt[2] - legend.width
vmar = (legend.plt[4] - legend.plt[3]) * ((1 - legend.params$shrink)/2)
legend.plt[4] = legend.plt[4] - vmar
legend.plt[3] = legend.plt[3] + vmar
image.plt[2] = min(image.plt[2], legend.plt[1] - hoffset)
par(plt = legend.plt)
key.marks = round(seq(from = 0, to = 1, length.out = legend.params$n.ticks) *
n.levels)
key.labels = formatC(as.numeric(wavelet.levels), digits = legend.params$label.digits,
format = legend.params$label.format)[key.marks +
1]
image(1, seq(from = 0, to = n.levels), matrix(wavelet.levels,
nrow = 1), col = key.cols, breaks = wavelet.levels,
useRaster = T, xaxt = "n", yaxt = "n", xlab = "",
ylab = "")
axis(4, lwd = lwd.axis, at = key.marks, labels = NA,
tck = 0.02, tcl = (par()$usr[2] - par()$usr[1]) *
legend.params$width - 0.04)
mtext(key.labels, side = 4, at = key.marks, line = 0.5,
las = 2, font = par()$font.axis, cex = par()$cex.axis)
text(x = par()$usr[2] + (1.5 + legend.params$lab.line) *
par()$cxy[1], y = n.levels/2, labels = legend.params$lab,
xpd = NA, srt = 270, font = par()$font.lab, cex = par()$cex.lab)
box(lwd = lwd.axis)
par(new = TRUE, plt = image.plt)
}
image(axis.1, axis.2, t(Power), col = key.cols, breaks = wavelet.levels,
useRaster = useRaster, ylab = "", xlab = "", axes = FALSE,
main = main)
if ((plot.contour == T) & (is.null(Power.pval) == F)) {
contour(axis.1, axis.2, t(Power.pval) < siglvl, levels = 1,
lwd = lwd, add = TRUE, col = col.contour, drawlabels = FALSE)
}
if (plot.ridge == T) {
Ridge = Ridge * (Power >= lvl)
contour(axis.1, axis.2, t(Ridge), levels = 1, lwd = lwd,
add = TRUE, col = col.ridge, drawlabels = FALSE)
}
if (plot.coi == T) {
polygon(WT$coi.1, WT$coi.2, border = NA, col = rgb(1,
1, 1, 0.5))
}
box(lwd = lwd.axis)
if (label.period.axis == T) {
period.axis.default = (is.null(spec.period.axis$at))
if (!period.axis.default) {
if (is.numeric(spec.period.axis$at)) {
period.axis.warning = ((sum(spec.period.axis$at <=
0, na.rm = T) > 0) | (sum(!is.na(spec.period.axis$at)) !=
sum(is.finite(spec.period.axis$at))))
}
else {
period.axis.warning = TRUE
}
if (is.logical(spec.period.axis$labels)) {
period.axis.warning = (period.axis.warning |
ifelse(length(spec.period.axis$labels) != 1,
TRUE, is.na(spec.period.axis$labels)))
}
if (!is.logical(spec.period.axis$labels)) {
period.axis.warning = (period.axis.warning |
(length(spec.period.axis$labels) != length(spec.period.axis$at)))
}
}
period.axis.default = (period.axis.default | period.axis.warning)
if ((is.null(periodlab)) | (!is.null(periodlab) & period.axis.warning)) {
periodlab = "period"
}
if (period.axis.default) {
period.tick = unique(trunc(axis.2))
period.tick[period.tick < log2(WT$Period[1])] = NA
period.tick = na.omit(period.tick)
period.tick.label = 2^(period.tick)
axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
tck = periodtck, tcl = periodtcl)
axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
tck = periodtck, tcl = periodtcl)
mtext(period.tick.label, side = 2, at = period.tick,
las = 1, line = par()$mgp[2] - 0.5, font = par()$font.axis,
cex = par()$cex.axis)
}
if (!period.axis.default) {
period.tick = log2(spec.period.axis$at)
period.tick[(period.tick < log2(WT$Period[1]))] = NA
period.tick.label = spec.period.axis$labels
which.na = which(is.na(period.tick))
if (length(which.na) > 0) {
period.axis.warning.na = T
}
if (is.logical(period.tick.label)) {
if (period.tick.label == T) {
period.tick.label = 2^(period.tick)
}
}
axis(2, lwd = lwd.axis, at = period.tick, labels = period.tick.label,
tck = periodtck, tcl = periodtcl, las = spec.period.axis$las,
hadj = spec.period.axis$hadj, padj = spec.period.axis$padj,
mgp = par()$mgp - c(0, 0.5, 0), font = par()$font.axis,
cex.axis = par()$cex.axis)
axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
tck = periodtck, tcl = periodtcl)
}
mtext(periodlab, side = 2, line = par()$mgp[1] - 0.5,
font = par()$font.lab, cex = par()$cex.lab)
}
if (label.time.axis == T) {
time.axis.default = (is.null(spec.time.axis$at))
if (show.date == T) {
if (is.element("date", names(series.data))) {
my.date = series.data$date
}
else {
my.date = rownames(series.data)
}
if (is.null(date.format)) {
chronology.warning = inherits(try(as.Date(my.date,
tz = date.tz), silent = T), "try-error")
if (!chronology.warning) {
my.date = as.Date(my.date, tz = date.tz)
chronology.warning = ifelse(sum(is.na(my.date)) >
0, TRUE, sum(diff(my.date, tz = date.tz) <
0) > 0)
}
}
if (!is.null(date.format)) {
chronology.warning = inherits(try(as.POSIXct(my.date,
format = date.format, tz = date.tz), silent = T),
"try-error")
if (!chronology.warning) {
my.date = as.POSIXct(my.date, format = date.format,
tz = date.tz)
chronology.warning = ifelse(sum(is.na(my.date)) >
0, TRUE, sum(diff(my.date, tz = date.tz) <
0) > 0)
}
}
if (chronology.warning) {
show.date = F
time.axis.default = T
timelab = "index"
}
}
if (!time.axis.default) {
if (show.date == F) {
time.axis.warning = (!is.numeric(spec.time.axis$at))
}
if (show.date == T) {
if (is.null(date.format)) {
time.axis.warning = inherits(try(as.Date(spec.time.axis$at,
tz = date.tz), silent = T), "try-error")
if (!time.axis.warning) {
time.axis.warning = (sum(!is.na(as.Date(spec.time.axis$at,
tz = date.tz))) == 0)
}
}
if (!is.null(date.format)) {
time.axis.warning = inherits(try(as.POSIXct(spec.time.axis$at,
format = date.format, tz = date.tz), silent = T),
"try-error")
if (!time.axis.warning) {
time.axis.warning = (sum(!is.na(as.POSIXct(spec.time.axis$at,
format = date.format, tz = date.tz))) ==
0)
}
}
}
if (!time.axis.warning) {
if (is.logical(spec.time.axis$labels)) {
time.axis.warning = (ifelse(length(spec.time.axis$labels) !=
1, TRUE, is.na(spec.time.axis$labels)))
}
if (!is.logical(spec.time.axis$labels)) {
time.axis.warning = (length(spec.time.axis$labels) !=
length(spec.time.axis$at))
}
}
}
time.axis.default = (time.axis.default | time.axis.warning)
}
half.increment.2 = (axis.2[2] - axis.2[1])/2
my.ylim = c(min(axis.2) - half.increment.2, max(axis.2) +
half.increment.2)
index.condition = ((show.date == F) | (label.time.axis ==
F))
if ((WT$dt != 1) & (index.condition)) {
par(new = TRUE)
plot(1:WT$nc, seq(min(axis.2), max(axis.2), length.out = WT$nc),
xlim = c(0.5, WT$nc + 0.5), ylim = my.ylim, type = "n",
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", xlab = "",
ylab = "")
}
if (label.time.axis == T) {
if ((is.null(timelab) & time.axis.default) | (!is.null(timelab) &
time.axis.warning)) {
timelab = ifelse(show.date, "calendar date", "index")
}
if (show.date == F) {
if (time.axis.default) {
A.1 = axis(1, lwd = lwd.axis, labels = NA, tck = timetck,
tcl = timetcl)
mtext(A.1, side = 1, at = A.1, line = par()$mgp[2] -
0.5, font = par()$font.axis, cex = par()$cex.axis)
}
if (!time.axis.default) {
time.tick = spec.time.axis$at
time.tick.label = spec.time.axis$labels
which.na = which(is.na(time.tick))
if (length(which.na) > 0) {
time.axis.warning.na = T
}
axis(1, lwd = lwd.axis, at = time.tick, labels = time.tick.label,
tck = timetck, tcl = timetcl, las = spec.time.axis$las,
hadj = spec.time.axis$hadj, padj = spec.time.axis$padj,
mgp = par()$mgp - c(0, 0.5, 0), font = par()$font.axis,
cex.axis = par()$cex.axis)
}
mtext(timelab, side = 1, line = par()$mgp[1] - 1,
font = par()$font.lab, cex = par()$cex.lab)
}
if (show.date == T) {
half.increment = difftime(my.date[2], my.date[1],
tz = date.tz)/2
if (is.null(date.format)) {
my.xlim = as.Date(range(my.date) - c(half.increment,
-half.increment), tz = date.tz)
}
else {
my.xlim = as.POSIXct(range(my.date) - c(half.increment,
-half.increment), format = date.format, tz = date.tz)
}
par(new = TRUE)
if (time.axis.default) {
plot(my.date, seq(min(axis.2), max(axis.2), length.out = WT$nc),
xlim = my.xlim, ylim = my.ylim, type = "n",
xaxs = "i", yaxs = "i", yaxt = "n", xlab = "",
ylab = "", lwd = lwd.axis, tck = timetck, tcl = timetcl,
mgp = par()$mgp - c(0, 0.5, 0), font = par()$font.axis,
cex.axis = par()$cex.axis)
}
if (!time.axis.default) {
plot(my.date, seq(min(axis.2), max(axis.2), length.out = WT$nc),
xlim = my.xlim, ylim = my.ylim, type = "n",
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")
if (is.null(date.format)) {
time.tick = as.Date(spec.time.axis$at, tz = date.tz)
}
if (!is.null(date.format)) {
time.tick = as.POSIXct(spec.time.axis$at, format = date.format,
tz = date.tz)
}
time.tick.label = spec.time.axis$labels
which.na = which(is.na(time.tick))
if (length(which.na) > 0) {
time.axis.warning.na = T
}
if (is.logical(time.tick.label)) {
if (time.tick.label == T) {
time.tick.label = time.tick
}
}
axis(1, lwd = lwd.axis, at = time.tick, labels = time.tick.label,
tck = timetck, tcl = timetcl, las = spec.time.axis$las,
hadj = spec.time.axis$hadj, padj = spec.time.axis$padj,
mgp = par()$mgp - c(0, 0.5, 0), font = par()$font.axis,
cex.axis = par()$cex.axis)
}
mtext(timelab, side = 1, line = par()$mgp[1] - 1,
font = par()$font.lab, cex = par()$cex.lab)
}
}
options(default.options)
if (graphics.reset == T) {
par(op)
}
if (period.axis.warning == T) {
warning("\nPlease check your period axis specifications. Default settings were used.")
}
if (period.axis.warning.na == T) {
warning("\nNAs were produced with your period axis specifications.")
}
if (chronology.warning) {
warning("\nPlease check your calendar dates, format and time zone: dates may not be in an unambiguous format or chronological. The default numerical axis was used instead.")
}
if (time.axis.warning == T) {
warning("\nPlease check your time axis specifications. Default settings were used.")
}
if (time.axis.warning.na == T) {
warning("\nNAs were produced with your time axis specifications.")
}
output = list(op = op, image.plt = image.plt, legend.plt = legend.plt)
class(output) = "graphical parameters"
out("Class attributes are accessible through following names:\n")
out(names(output), "\n")
return(invisible(output))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.