Train yourself by reproducing this code in python. For in-depth analyses, I strongly recommend the use of the ipython notebook. This tutorial was made with that tool.
Start importing wavepal and other usual packages
In [1]:
% matplotlib inline
# other option (instead of "inline"): nbagg
# N.B.: the two previous lines are used by ipython. Do not write them in your python code.
import numpy as np
import matplotlib.pyplot as plt
import copy
import wavepal as wv
wavepal is illustrated here with ODP1148 d18O benthic time series. Ref: Z. Jian, Q. Zhao, X. Cheng, J. Wang, P. Wang, and X. Su. Pliocene-pleistocene stable isotope and paleoceanographic changes in the northern south china sea. Palaeogeography, Palaeoclimatology, Palaeoecology, 193(3–4):425–442, 2003.
Read the data
In [2]:
data=np.genfromtxt("ODP1148-BF-18O.txt")
myt=data[:,0]
mydata=data[:,1]
Initialize the class called Wavepal (which is a class of the package wv)
In [3]:
x=wv.Wavepal(myt, mydata, "Age", "$\delta{}^{18}O$", t_units="ka", mydata_units="permil")
To get more info about the inputs at the initialization, use the help:
In [4]:
help(x.__init__)
Here is the list of all the methods of the Wavepal class:
For each method, you have at your disposal a detailed explanation about what it does and about the inputs and outputs, thanks to "help". Example with freq_analysis method:
In [5]:
help(x.freq_analysis)
N.B.: you can get all of them with: help(x.__module__)
Below is presented an example of spectral analysis using some of those methods.
Check the data set
In [6]:
x.check_data()
Figure of the time step in function of time, with an histogram of the distribution of the time steps. If you want to save the figure, use "savefig".
In [7]:
plot_timestep=x.plot_timestep()
plot_timestep.savefig("timestep.pdf")
plot_timestep.show()
Figure of the trend. Try many degrees of the polynomial.
In [8]:
plot_trend=x.plot_trend(pol_degree=-1) # no trend
plot_trend.show()
plot_trend=x.plot_trend(pol_degree=0) # constant trend
plot_trend.show()
plot_trend=x.plot_trend(pol_degree=5)
plot_trend.show()
plot_trend=x.plot_trend(pol_degree=7)
plot_trend.show()
plot_trend=x.plot_trend(pol_degree=9)
plot_trend.show()
and choose the degree of the polynomial for the subsequent analyses
In [9]:
x.choose_trend_degree(7)
Compute some variables related to the trend
In [10]:
x.trend_vectors()
Option make_carma_fig=True allows to save some figures related to the analyses of the background noise (made with 'carma pack' package of Kelly & al. - see the reference in the help of carma_params method). They are here saved in the folder "figures_carma", previously created by the user.
In [11]:
x.carma_params(make_carma_fig=True,nbins=20,dpi=400,path_to_figure_folder="figures_carma/")
Choose the $\textrm{x}^{\textrm{th}}$ percentiles for significance testing
In [12]:
percentile=np.zeros(2)
percentile[0]=95.
percentile[1]=99.9
N.B.: If you are not interested in the frequency analysis, you can skip the code below and go to section 5.
Run the method for the frequency analysis (do not forget to use help(x.freq_analysis) if needed)
In [13]:
x.freq_analysis(freqstep=0.0001,D=600.,percentile=percentile,n_moments=12,computes_amplitude=True)
Periodogram. This is the central result for the frequency analysis.
In [14]:
plot_periodogram=x.plot_periodogram(fontsize_legend=10)
plot_periodogram.show()
Figure of the WOSA segmentation
In [15]:
plot_WOSA_segmentation=x.plot_WOSA_segmentation()
plot_WOSA_segmentation.show()
Figure of the number of WOSA segments for each frequency
In [16]:
plot_number_WOSA_segments=x.plot_number_WOSA_segments()
plot_number_WOSA_segments.show()
Check the convergence of the analytical confidence levels
In [17]:
plot_check_convergence_percentiles=x.plot_check_convergence_percentiles(fontsize_suptitle=13,fontsize_title=9,fontsize_axes=8,fontsize_ticks=6)
plot_check_convergence_percentiles.show()
Compare the squared amplitude vs the periodogram
In [18]:
plot_amplitude_vs_periodogram=x.plot_amplitude_vs_periodogram(fontsize_legend=10)
plot_amplitude_vs_periodogram.show()
Define the times at which the CWT (continuous wavelet transform) is to be computed
In [19]:
theta=np.linspace(myt[0],myt[-1],2000)
Run the method for the time-frequency analysis. Note that the percentages of the progressing bars are not fully reliable, since the iterations in those loops do not spend equal computing time.
In [20]:
x.timefreq_analysis(theta=theta,w0=5.5,permin=10.,percentile=percentile)
Set the location of the ticks
In [21]:
time_string=[0., 500., 1000., 1500., 2000., 2500., 3000., 3500., 4000., 4500., 5000., 5500., 6000.]
period_string=[10., 21., 41., 100., 200., 400., 800., 1500.]
dashed_periods=[21., 41., 100.]
Scalogram. This is the central result for the time-frequency analysis. Resizing (with set_size_inches) is better, and is compulsory before saving the figure (with savefig).
In [22]:
plot_scalogram=x.plot_scalogram(color_cl_anal=['g'],color_cl_mcmc=['m'],fontsize_title=30,fontsize_axes=20,fontsize_ticks=20,linewidth_cl=4,global_scal_xlabel_ticks="bottom",decimals=2,linewidth_gscal=2.0,time_string=time_string,period_string=period_string,dashed_periods=dashed_periods,power_string=power_string)
plot_scalogram=x.plot_scalogram(color_cl_anal=['m','g'],fontsize_title=40,fontsize_ticks=20,fontsize_axes=30,time_string=time_string,period_string=period_string,dashed_periods=dashed_periods)
fig = plt.gcf()
fig.set_size_inches(24,12)
plot_scalogram.show()
Same analysis, with another time-frequency resolution (change the "w0" parameter of the Morlet Wavelet)
In [23]:
y=copy.copy(x) # optional: in order to keep in memory all what we did up to here (very convenient if you work with ipython)
y.timefreq_analysis(theta=theta,w0=15.0,permin=10.,percentile=percentile)
In [24]:
plot_scalogram=y.plot_scalogram(color_cl_anal=['m','g'],fontsize_title=40,fontsize_ticks=20,fontsize_axes=30,time_string=time_string,period_string=period_string,dashed_periods=dashed_periods)
fig = plt.gcf()
fig.set_size_inches(24,12)
plot_scalogram.show()
You may want to access the attributes of the Wavepal class. The list of available attributes is given in the method __init__ of the Wavepal class (open the file Wavepal.py). Most of them are initialized as None. Their values may change when running the methods of Wavepal. Consequently, I recommend to read, save, and possibly modify (with caution) for further use, the attributes of Wavepal after you have run all the desired methods.
Getting the value of an attribute is easy. Example with 'nt' (number of data points):
In [25]:
x.nt
Out[25]:
The attributes are listed below, with some details given.
Two methods that are not part of the Wavepal class are also available to the user, in the wv package:
In [26]:
help(wv.mov_av)
In [27]:
help(wv.dt_GCD)