The following notebook contains examples for using the psprocess.py and then the pscrosscorr toolbox for preprocessing raw seismic waveforms. These are then cross-correlated to produce a Green's function, representative of the trace path between the two stations. Currently the script can only operate with MSEED formats, but additional support for other formats such as: SAC, SEED and SUDS will likely be added in the future. The user specifies the input path to the raw waveform. This waveform should be only one trace for the purposes of this example, e.g. BHZ. A Preprocess object is created with this input, and the output is specified by one of the functions contained.
The theory for the current pre-processing workflow in this example follows for preprocessing is explained in depth in Bensen et al. (2007). This example is a simple one that shows the outputs of two cross-correlated, unstacked, pre-processed seismic waveforms. If you want to learn more about the pre-processing steps, please refer to the ipython notebook examples provided.
In [8]:
from pysismo.pspreprocess import Preprocess
from obspy.signal.cross_correlation import xcorr
from obspy import read
from obspy.core import Stream
import matplotlib.pyplot as plt
import numpy as np
import os
%matplotlib inline
The Preprocess class requires many input parameters to function. Below is a list of examples.
In [12]:
# list of example variables for Preprocess class
FREQMAX = 1./1 # bandpass parameters
FREQMIN = 1/20.0
CORNERS = 2
ZEROPHASE = True
ONEBIT_NORM = False # one-bit normalization
PERIOD_RESAMPLE = 0.02 # resample period to decimate traces, after band-pass
FREQMIN_EARTHQUAKE = 1/75.0 # earthquakes periods band
FREQMAX_EARTHQUAKE = 1/25.0
WINDOW_TIME = 0.5 * 1./FREQMAX_EARTHQUAKE # time window to calculate time-normalisation weights
WINDOW_FREQ = 0.0002 # freq window (Hz) to smooth ampl spectrum
CROSSCORR_TMAX = 300 # set maximun cross-correlation time window both positive and negative!
# initialise the Preprocess class
PREPROCESS = Preprocess(FREQMIN, FREQMAX, FREQMIN_EARTHQUAKE,
FREQMAX_EARTHQUAKE, CORNERS, ZEROPHASE,
PERIOD_RESAMPLE, WINDOW_TIME, WINDOW_FREQ,
ONEBIT_NORM)
In [7]:
# specify the path to the desired waveforms, the example folder xcorr is provided.
example_folder = 'tools/examples/xcorr'
def paths(folder_path, extension, sort=False):
"""
Function that returns a list of desired absolute paths called abs_paths
of files that contains a given extension e.g. .txt should be entered as
folder_path, txt. This function will run recursively through and find
any and all files within this folder with that extension!
"""
abs_paths = []
for root, dirs, files in os.walk(folder_path):
for f in files:
fullpath = os.path.join(root, f)
if os.path.splitext(fullpath)[1] == '.{}'.format(extension):
abs_paths.append(fullpath)
if sort:
abs_paths = sorted(abs_paths, key=paths_sort)
return abs_paths
# create a list of possible MSEED files inside this folder to iterate over.
abs_paths = paths(example_folder, 'mseed')
preprocessed_traces = []
for abs_path in abs_paths:
st_headers = read(abs_path, headonly=True)
starttime = st_headers[0].stats.starttime
endtime = starttime + 10*60 #only input 10min waveforms for this example
# import a trace from the example waveform
st = read(abs_path, starttime=starttime, endtime=endtime)
# merge stream
st.merge(fill_value='interpolate')
# select trace from merged stream
tr = st[0]
#preprocess the trace
tr = PREPROCESS.bandpass_filt(tr)
tr = PREPROCESS.trace_downsample(tr)
# copy the trace for time normalisation
tr_copy = tr
# process for time normalisation
tr = PREPROCESS.time_norm(tr, tr_copy)
tr = PREPROCESS.spectral_whitening(tr)
preprocessed_traces.append(tr)
Now perform a cross-correlation using obspy's c-wrapped xcorr function.
In [17]:
# set the two traces tr1 and tr2 to be cross-correlated
tr1, tr2 = preprocessed_traces
# initializing time and data arrays of cross-correlation, because xcorrs can be symmetric, the final waveform will be
# symmetic about 0 where the maximum extents are CROSSCORR_TMAX in seconds both +ve and -ve.
nmax = int(CROSSCORR_TMAX / PERIOD_RESAMPLE)
timearray = np.arange(-nmax * PERIOD_RESAMPLE, (nmax + 1)*PERIOD_RESAMPLE, PERIOD_RESAMPLE)
# the shift length provides the ability to choose the overlap between the traces to be correlated.
SHIFT_LEN = (len(timearray) - 1) * 0.5
crosscorr = xcorr(tr1, tr2, shift_len=len(tr1.data)/2, full_xcorr=True)[2]
plt.figure()
plt.plot(crosscorr)
plt.show()