Another test on my MacBook Pro: The wavelet transform example in Megies et al (2011).
You may need to install the Machine Learning Python package. The command
conda install mlpy
did not work for me, so I search the internet for a method to install
mlpy
The mlpy package requires the GNU Scientific Library (GSL), which was not installed on my laptop, and may not be installed on yours either. This package can be found on the web. My method was taken from the book "Practical Data Analysis Cookbook" by Tomasz Drabas. First, direct your web browser to
http://gnu.mirro.iweb.com/gsl
and download the latest version of GSL, typically called
gsl-latest.tar-gz. Currently the latest is V2.4. Move the downloaded file into your Applications directory (on a Mac) -- I dont know the equivalent place on a Windows machine. Then unpack this gzipped tarfile by double-clicking on it ot else use the command
gzip -c gsl-latest.tar-gz | tar -xv
The unpacked file will produce a folder with the GSL software in it. In your terminal-app go to this folder and prepare and install GSL
./configure
make
make install
Hopefully this compilation and installation proceeds without error. Now download a compressed tarfile with the lastest version of
mlpy from this website into your
anaconda dirctory, which on my laptop is
/Users/jjpark/anaconda. The current latest version is V3.5.
http://sourceforge.net/projects/mlpy/files
Double-click the tarfile to extract its contents into a subdirectory, or else
gzip -c mlpy-.tar.gz | tar -xv
or, if already uncompressed,
tar -xvf mlpy-.tar
This action produces a subdirectory
mlpy-.
cd to this directory in your terminal app and type
python setup.py install
and
mlpy will be installed if all goes well. Python modules that one obtains from sourceforge and other websites typically include a
setup.py scriptfile to instruct python how to set the module up for later use.
I discovered coding bugs in the
mlpy module while trying to execute the coding examples. The continuous wavelet transform relies on a script file
./mlpy-3.5.0/mlpy/wavelet/continuous.py that had a mistake in its Line 162:
161
J = floor(dj**-1 * log2((N * dt) / s0))
162
s = empty(J + 1)
The python function "
floor" from the
numpy module returns a floating-point number, but the python function "
empty" expects an integer to define how long an array to generate. This combination of mismatched variable types leads to an error message:
'numpy.float64' object cannot be interpreted as an integer
Fixing the code to read, e.g.,
J = floor(dj**-1 * log2((N * dt) / s0))
JJ = int(J)
s = empty(JJ + 1)
*should* solve the problem, but it doesn't. Why? Because the original
continuous.py is loaded into the python libraries during the
install-operation, and changing the code after installation doesnt change the installed code. You must make changes to the script
continuous.py in the
mlpy-3.5.0 subdirectory, and then re-install the module by executing
python setup.py install
while the terminal app is "in" the
./mlpy-3.5.0 directory. Note that python will tell you that the execution error resides in a copy of
continuous.py that lies deep within a subdirectory tree. This is NOT the copy that needs to be debugged. You must debug the copy of
continuous.py that sits close within the subdirectory
anaconda/mlpy-3.5.0/mlpy/wavelet, and then re-install the
mlpy module.
Here is a simple example of a wavelet transform of a randomly-generated array of numbers
import numpy as np
import matplotlib.pyplot as plt
import mlpy.wavelet as wave
x = np.random.sample(512)
scales = wave.autoscales(N=x.shape[0], dt=1, dj=0.025, wf='morlet', p=2)
X = wave.cwt(x=x, dt=1, scales=scales, wf='dog', p=2)
fig = plt.figure(1)
ax1 = plt.subplot(2,1,1)
p1 = ax1.plot(x)
ax1.autoscale_view(tight=True)
ax2 = plt.subplot(2,1,2)
p2 = ax2.imshow(np.abs(X), interpolation='nearest')
plt.show()
Here is the wavelet transform of the data series we downloaded in the first ObsPy exercise
import numpy as np, matplotlib.pyplot as plt
import mlpy
import mlpy.wavelet as wave
from obspy.core import UTCDateTime
from obspy.clients.arclink import Client
from obspy.signal.filter import lowpass
from obspy.signal.invsim import corn_freq_2_paz, simulate_seismometer
client = Client("webdc.eu", port=18001)
t = UTCDateTime("2009-08-24 00:20:03")
st = client.get_waveforms("BW", "RJOB", "", "EHZ", t, t+30)
omega0 = 2.0
trace = st[0].data
npts = len(trace)
dt = st[0].stats.delta
dj = 0.05
wf = 'morlet'
p=omega0
scl = mlpy.wavelet.autoscales(npts, dt, dj , wf, p)
spec = mlpy.wavelet.cwt(trace, dt, scl , wf='morlet', p=omega0)
freq = (omega0 + np.sqrt(2.0 + omega0**2)) / (4*np.pi * scl[1:])
t = np.arange(st[0].stats.npts) / st[0].stats.sampling_rate
plt.imshow(np.abs(spec), extent=(t[0], t[-1], freq[-1], freq[0]))
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.show()
Try these scripts to be sure that your installation of Python and ObsPy is working properly.
