wiki:RoseMilan
Last modified 16 months ago

Work in progress

1 Intro ObsPy/Python

1.1 Basic Matrix Computation

This example shows some basic matrix computation using the numpy.ndarray object ( http://numpy.scipy.org )
. Example Script:  math_python.py  Download

Exercise:

  • Calculate maximum, minimum, mean, std of the array x. Find available methods on the numpy.ndarray here:  http://numpy.scipy.org.

2 Data formats

3 ArcLink, DHI/Fissures

Retrieve seismograms via WebDC with the ArcLink protocol.
Example Script:  getwaveform_arclink.py  Download

Exercise:

  • Retrieve data from different stations. Stations and network ids can be found here:  http://www.webdc.eu -> Query
  • Try different channels (also available at  http://www.webdc.eu), determine sampling_rate of data from different channel
  • Directly save waveform (uncomment last part in script), create data directory and move saved waveforms in it

4 Instrument Correction

The response of a seismometer can be written as poles and zeros in the Laplace domain. In the following example the instrument response of a Lenarz LE3D seismometer is deconvolved and a Wood Anderson seismometer is simulated.
Example Script:  correct_instrument.py  Download
 Download Data

Exercise:

  • XXX: Heiner Magnitude
  • Correct only the instrument without simulating an other one. Therefore us inst_sim=None. This will result in displacement [nm]
  • Correct to [nm/s] by integrating the signal. This is accomplished by removing one zero in the le3d PAZ.

5 Spectral Estimation

5.1 Tapering

In this example we estimate the frequency of a sine wave. As the sine wave is cut in such a way that it is not mirror side repeatable, there is loss to other frequencies, a taper can circumvent this.
Example Script:  tapering.py  Download

Exercise:

  • Change cut off frequency of bandpass.
  • Change generator frequency of sine.
  • Use blackman taper window

5.2 Eigenvalues

5.3 Downsampling

As data two sinus are taken with 4Hz and 8Hz. In order to downsample the data, every second item is taken once with zeroing frequencies above the nyquist frequency an once without.
Example Script:  downsampling.py  Download

Exercise:

  • Change the generator frequencies
  • Instead of zeroing frequencies above Nyquist lowpass filter the data with Nyquist frequency

5.4 Detrending

5.5 Filtering

5.5.1 Filtering with a simple boxcar

A simple filter is constructed by transforming the signal to the frequency domain and zeroing all unwanted frequencies. Study the construction of a filter given in the function  filter_bc.py
Example Scripts:  filter_bc.py  Download;  simple_filters.py  Download

Exercise:

  • Try different cut off frequencies
  • Take filter_bc.py as a guideline and construct a simple bandpass in a similar way

5.5.2 Filtering with an Butterworth filter

A standard lowpass filter is used to filter a artificial spike and a real signal:
Example Scripts:  filter_bw_lp.py  Download  Download germany.mat

Exercise:

  • See effects for different cut off frequencies f0
  • Copy the script and do the same with highpass and bandpass filters

5.6 Integration to Displacement

5.7 Generate synthetic Seismogram

Generate a seismogram by convolving a Greens function with a source time function. Spikes at (-1,1) are taken as representation of the Greens function, the source time function is represented as a ricker wavelet.
Example Script:  generate_seismogram.py  Download

Exercise:

  • Play with the parameters of the ricker wavelet.
  • Set more/less spikes in the Greens function.









Effect of tapers for spectrum (arbitrary amplitude)

Import necessary modules: 

from obspy.core import read
from obspy import UTCDateTime
from copy import deepcopy
import matplotlib
from matplotlib.pyplot import *
from numpy import *

Read in the data, the data can be downloaded here. 

tr_orig = read("loc_RNON20040609200559.z")[0]
tr = deepcopy(tr_orig)
figure(1)
time = linspace(0, tr.stats.npts/tr.stats.sampling_rate, tr.stats.npts)
plot(time, tr.data, label='full length data', linewidth=1)
legend()
xlabel('Time [s]')

Zoom the data to appropriate length 

T1 = UTCDateTime("20040609200621.12")
T2 = UTCDateTime("20040609200623.79")
tr.trim(T1,T2)
npts, nyq = tr.stats.npts, tr.stats.sampling_rate/2.0 #nsamp, nyquist
time = linspace(0, npts/nyq/2, npts)
figure(2)
plot(time, tr.data, label='trimmed trace')
xlabel('Time [s]')

Fourier transform 

tr_f = fft.rfft(tr.data) #real fft
freq = linspace(0,nyq,npts/2)
figure(3)
plot(freq, abs(tr_f[1:]), label='no taper')

Different tapers 

from obspy.signal import cosTaper
win_ham = hamming(npts)
win_cos = cosTaper(npts,0.1) #10%
figure(4)
plot(time, win_ham, label='hamming')
plot(time, win_cos, label='10% cosine')
ylim(0,2)
legend()
xlabel('Time [s]')

Taper the signal 

figure(2)
plot(time, tr.data * win_ham, label='hamming')
plot(time, tr.data * win_cos, label='10% cosine')
legend()
xlabel('Time [s]')

Spectrum of tapered signal (arbitrary y scale) 

tr_fh = fft.rfft(tr.data*win_ham)
tr_fc = fft.rfft(tr.data*win_cos)
figure(3)
plot(freq, abs(tr_fh[1:]), label='hamming')
plot(freq, abs(tr_fc[1:]), label='10% cosine')
ax = gca()
#ax.set_xscale('log')
ax.set_yscale('log')
legend()
xlabel('Frequency [Hz]')
show()

Built in power spectrum which also includes detrending and normalization due to windowing loss. 

from matplotlib import psd
figure()
psd(tr.data,npts+1,Fs=nyq*2)

Attachments