This document describes how to perform a wavelet analysis of an annual precipitation timeseries. Three identical analyses are performed using 1) Scott Steinshneider's code, and 2) the biwavelet
package, and 3) the weathergen
package.
library(dplyr) library(tidyr) library(lubridate) library(ggplot2) theme_set(theme_bw()) library(weathergen) data(climate) # compute daily averages across sites clim.da <- group_by(climate, DATE) %>% summarise(N=n(), PRCP=mean(PRCP), TMIN=mean(TMIN), TMAX=mean(TMAX)) # aggregate by water year clim.wyr <- clim.da %>% mutate(WYEAR=wyear(DATE)) %>% group_by(WYEAR) %>% summarise(N=n(), PRCP=sum(PRCP), TMIN=mean(TMIN), TMAX=mean(TMAX))
Here is a timeseries plot of the annual precipitation timeseries.
ggplot(clim.wyr, aes(WYEAR, PRCP)) + geom_line() + labs(x="Water Year", y="Annual Precipitation (mm)")
The goal of the wavelet analysis is to identify any significant low-frequency signals in the annual precipitation timeseries.
This section uses the code written by Scott Steinschneider.
WAVELET_ANALYSIS <- function(CLIMATE_VARIABLE,siglvl,background_noise,plot_flag) { if (missing(plot_flag)) {plot_flag <- TRUE} if (missing(siglvl)) {siglvl <- 0.90} if (missing(background_noise)) {background_noise <- "white"} #Define Wavelet function used that returns transform waveletf <- function(k,s) { nn <- length(k) k0 <- 6 #nondimensional frequency, here taken to be 6 to satisfy the admissibility condition [Farge 1992] z <- array(1,nn) z[which(k<=0)] <- 0 expnt <- -((s*k - k0)^2/2)*z norm <- sqrt(s*k[2])*(pi^(-0.25))*sqrt(nn) # total energy=N [Eqn(7)] daughter <- norm*exp(expnt) daughter <- daughter*z fourier_factor <- (4*pi)/(k0 + sqrt(2 + k0^2)) # Scale-->Fourier [Sec.3h] coi <- fourier_factor/sqrt(2) # Cone-of-influence [Sec.3g] dofmin <- 2 # Degrees of freedom return(daughter) } #Define Wavelet function used that returns fourier_factor, cone of influence, and degrees of freedom waveletf2 <- function(k,s) { nn <- length(k) k0 <- 6 #nondimensional frequency, here taken to be 6 to satisfy the admissibility condition [Farge 1992] z <- array(1,nn) z[which(k<=0)] <- 0 expnt <- -((s*k - k0)^2/2)*z norm <- sqrt(s*k[2])*(pi^(-0.25))*sqrt(nn) # total energy=N [Eqn(7)] daughter <- norm*exp(expnt) daughter <- daughter*z fourier_factor <- (4*pi)/(k0 + sqrt(2 + k0^2)) # Scale-->Fourier [Sec.3h] coi <- fourier_factor/sqrt(2) # Cone-of-influence [Sec.3g] dofmin <- 2 # Degrees of freedom return(c(fourier_factor,coi,dofmin)) } # Perform Wavelet #....construct time series to analyze, pad if necessary CURRENT_CLIMATE_VARIABLE_org <- CLIMATE_VARIABLE variance1 <- var(CURRENT_CLIMATE_VARIABLE_org) n1 <- length(CURRENT_CLIMATE_VARIABLE_org) CURRENT_CLIMATE_VARIABLE <- scale(CURRENT_CLIMATE_VARIABLE_org) variance2 <- var(CURRENT_CLIMATE_VARIABLE) base2 <- floor(log(n1)/log(2) + 0.4999) # power of 2 nearest to N CURRENT_CLIMATE_VARIABLE <- c(CURRENT_CLIMATE_VARIABLE,rep(0,(2^(base2+1)-n1))) n <- length(CURRENT_CLIMATE_VARIABLE) #Determine parameters for Wavelet analysis dt <- 1 dj <- 0.25 s0 <- 2*dt J <- floor((1/dj)*log((n1*dt/s0),base=2)) #....construct SCALE array & empty PERIOD & WAVE arrays scale <- s0*2^((0:J)*dj) period <- scale wave <- array(as.complex(0),c(J+1,n)) # define the wavelet array #....construct wavenumber array used in transform [Eqn(5)] k <- c(1:floor(n/2)) k <- k*((2.*pi)/(n*dt)) k <- c(0,k,-rev(k[1:floor((n-1)/2)])) f <- fft(CURRENT_CLIMATE_VARIABLE,inverse=FALSE) #fourier transform of standardized precipitation # loop through all scales and compute transform for (a1 in 1:(J+1)) { daughter <- waveletf(k,scale[a1]) results <- waveletf2(k,scale[a1]) fourier_factor <- results[1] coi <- results[2] dofmin <- results[3] wave[a1,] <- fft(f*daughter,inverse=TRUE)/n # wavelet transform[Eqn(4)] } period <- fourier_factor*scale coi <- coi*dt*c((0.00001),1:((n1+1)/2-1),rev((1:(n1/2-1))),(0.00001)) # COI [Sec.3g] wave <- wave[,1:n1] # get rid of padding before returning POWER <- abs(wave)^2 GWS <- variance1*apply(POWER,FUN=mean,c(1)) #Global Wavelet Spectrum # Significance Testing # get the appropriate parameters [see Table(2)] k0 <- 6 empir <- c(2,0.776,2.32,0.60) dofmin <- empir[1] # Degrees of freedom with no smoothing Cdelta <- empir[2] # reconstruction factor gamma_fac <- empir[3] # time-decorrelation factor dj0 <- empir[4] # scale-decorrelation factor if (background_noise=="white") {lag1 <- 0} #for red noise background, lag1 autocorrelation = 0.72, for white noise background, lag1 autocorrelation = 0 if (background_noise=="red") {lag1 <- .72} freq <- dt / period # normalized frequency fft_theor <- (1-lag1^2) / (1-2*lag1*cos(freq*2*pi)+lag1^2) # [Eqn(16)] fft_theor <- fft_theor # include time-series variance dof <- dofmin #ENTIRE POWER SPECTRUM chisquare <- qchisq(siglvl,dof)/dof signif <- fft_theor*chisquare # [Eqn(18)] sig95 <- ((signif))%o%(array(1,n1)) # expand signif --> (J+1)x(N) array sig95 <- POWER / sig95 # where ratio > 1, power is significant #TIME_AVERAGED (GLOBAL WAVELET SPECTRUM) dof <- n1 - scale if (length(dof) == 1) {dof <- array(0,(J+1))+dof} dof[which(dof < 1)] <- 1 dof <- dofmin*sqrt(1 + (dof*dt/gamma_fac / scale)^2 ) tt <- which(dof < dofmin) dof[tt] <- dofmin chisquare_GWS <- array(NA,(J+1)) signif_GWS <- array(NA,(J+1)) for (a1 in 1:(J+1)) { chisquare_GWS[a1] <- qchisq(siglvl,dof[a1])/dof[a1] signif_GWS[a1] <- fft_theor[a1]*variance1*chisquare_GWS[a1] } # PLOT period_lower_limit <- 0 sig_periods <- which(GWS>signif_GWS & period > period_lower_limit) if (plot_flag) { par(mfrow=c(1,2),font.axis=2,font.lab=2,cex.axis=1.1,cex.lab=1.1) xx <- 1:n1 yy <- period image(x=xx,y=yy,z=t(log(POWER,base=2)),xlab="Time (Years)",ylab="Fourier Period (Years)",ylim=rev(range(yy)),log="y",col=(heat.colors(12))) lines(xx,coi,lty=2) contour(x=xx,y=yy,z=t(sig95),add=TRUE,levels=c(-99,1),labels="") xmin <- (min(GWS,signif_GWS)) xmax <- (max(GWS,signif_GWS)) plot(GWS,yy,type="b",xlim=c(xmin,xmax),xlab="Global Wavelet Spectrum",ylab="Fourier Period (Years)",log="y",ylim=rev(range(yy))) lines(signif_GWS,yy,col="red",lty=2) if (background_noise=="white") {tt <- paste(siglvl*100,"% confidence \nlevel for white-noise\n spectrum",sep="")} if (background_noise=="red") {tt <- paste(siglvl*100,"% confidence \nlevel for red-noise\n spectrum",sep="")} legend(.2*max(GWS),.65*max(yy),tt,col="red",lty=2,box.lty=0,box.lwd=0,cex=.8) mtext("SIGNIFICANT PERIODS:",line=2) e <- paste(sig_periods[1]) if (length(sig_periods)==0) { mtext("NONE!",line=1) } else if (length(sig_periods)==1) { mtext(e,line=1) } else { for (h in 2:length(sig_periods)) { e <- paste(e,",",paste(sig_periods[h]),sep="") } mtext(e,line=1) } } return(list(GWS,signif_GWS,period)) }
The wavelet analysis and significance test are performed using the WAVELET_ANALYSIS()
function. This function accepts a numeric array of the annual precipitation time series, a significance level, a type of background noise (white or red), and a flag for generating a plot of the result. It returns the global wave spectrum (GWS), significance levels, and associated period. The periods with powers greater than the corresponding significance levels are significant.
wave.scott <- WAVELET_ANALYSIS(CLIMATE_VARIABLE=clim.wyr$PRCP, siglvl=0.90, background_noise="white", plot_flag=TRUE)
The plots show 1) a heat map of the localized power (W(s,t)) and 2) the time-averaged global wavelet spectrum (GWS) as a function of Fourier period. The --o--o-- line shows fourier period vs. GWS power, and the red ----- line shows period vs. significance power. The "Significant Periods" are those where the GWS power > significance power.
The output of this function is an unnamed list of three vectors corresponding to the 1) GWS power, 2) significance levels, and 3) fourier periods:
wave.scott
If the significance level is lowered from 0.9 to 0.8 then the red line of significance shifts towards the y-axis.
WAVELET_ANALYSIS(CLIMATE_VARIABLE=clim.wyr$PRCP, siglvl=0.80, background_noise="white", plot_flag=TRUE)
The biwavelet
package contains two functions (biwavelet::wt()
and biwavelet::wt.sig()
) that provide identical functionality as Scott's WAVELET_ANALYSIS
function. The biwavelet
package is a translation of the original matlab code accompanying Torrence and Compo (1998). This same matlab code was used by Scott to create his own R version of the package. Therefore, both Scott's code and the biwavelet
package code are essentially identical.
Here is the same analysis for significance level of 0.90 and white background noise using the biwavelet
package.
library(biwavelet) bw <- biwavelet::wt(d=cbind(1:length(clim.wyr$PRCP), clim.wyr$PRCP), dt=1, dj=1/4, max.scale=length(clim.wyr$PRCP), sig.level=0.90, lag1=0, sig.test=1) str(bw)
This function returns an object of class biwavelet
that includes a generic plotting function.
par(mfrow=c(1,1)) plot(bw)
The time-averaged global wave spectrum (GWS) is then computed as the row-wise mean the $power
matrix.
bw$gws <- apply(bw$power, 1, mean)
The significance test of the GWS can be computed using the biwavelet::wt.sig()
function. Note that the degrees of freedom is equal to the number of data points adjusted by scale (to account for edge effects).
bw$gws.sig <- wt.sig(d=cbind(clim.wyr$WYEAR, clim.wyr$PRCP), dt=bw$dt, scale=bw$scale, sig.test=1, sig.level=0.90, dof=length(clim.wyr$PRCP)-bw$scale, mother='morlet', lag1=0.0)
The following tables show identical values for the period, global spectrum power, and significance between Scott's WAVELET_ANNUAL()
and the biwavelet::wt()
function output.
# compare periods cbind(biwavelet=bw$period, scott=wave.scott[[3]]) # compare global spectrum (row-wise mean of power for biwavelet::wt()) cbind(biwavelet=bw$gws, scott=wave.scott[[1]]) # compare signif cbind(biwavelet=bw$gws.sig$signif, scott=wave.scott[[2]])
We can also plot the output of biwavelet::wt()
par(mfrow=c(1,2)) plot(bw, plot.cb=TRUE, plot.phase=FALSE) plot(bw$gws, bw$period, type="b", xlab="Global Wavelet Spectrum", ylab="Fourier Period (Years)", log="y", ylim=rev(range(bw$period)), xlim=range(c(bw$gws, bw$gws.sig$signif))) lines(bw$gws.sig$signif, bw$period, lty=2, col='red')
Which is effectively identical to the WAVELET_ANALYSIS()
output plot.
WAVELET_ANALYSIS(CLIMATE_VARIABLE=clim.wyr$PRCP, siglvl=0.90, background_noise="white", plot_flag=TRUE)
The weathergen
package provides a wavelet_analysis()
function that uses the biwavelet wt()
and wt.sig()
functions to provide similar functionality as Scott's WAVELET_ANALYSIS()
function, but with fewer lines of code.
wave.wgen <- wavelet_analysis(x=clim.wyr$PRCP, years=clim.wyr$WYEAR, sig.level=0.90, noise.type="white") str(wave.wgen)
And we can plot this same as with the biwavelet
output as before.
par(mfrow=c(1,2)) plot(wave.wgen, plot.cb=TRUE, plot.phase=FALSE) plot(wave.wgen$gws, wave.wgen$period, type="b", xlab="Global Wavelet Spectrum", ylab="Fourier Period (Years)", log="y", ylim=rev(range(wave.wgen$period)), xlim=range(c(wave.wgen$gws, wave.wgen$gws.sig$signif))) lines(wave.wgen$gws.sig$signif, wave.wgen$period, lty=2, col='red')
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.