A collection of methods for data stacking

A collection of methods for stacking of generic time series data. See Yang et al. (2022)[https://academic.oup.com/gji/article/232/3/1600/6762921] for methods validation and comparison.

Please cite Yang et al. (2022)[https://academic.oup.com/gji/article/232/3/1600/6762921] if you used stackmaster functions in your research/publications.

1. Matlab: MATLAB implementations of the stacking methods. NOTE: tf-PWS in MATLAB is not ready yet and may produce unexpected stacking results.
2. stackmaster: root directory of the stackmaster codes in Python.
3. data: example data for test run.
4. figs: figures used in this README file.
5. notebooks: example notebooks to use the functions.

This package is under active development. The currently available modules are listed here.

1. utils: This module contains frequently used utility functions.

2. core: This module contains core stacking functions.

1. Install stackmaster package functions using pip. We recommend installing stackmaster in a virtual environment.

$ pip install stackmaster

This step will install the StackMaster package. The modules would then be imported under any working directory.

2. Install with local copy:

cd to the directory you want to save the package files. Then run:

$ pip install .

Run the following commands to test your installation, under the root directory of StackMaster.

iimport os,pickle
import numpy as np
import matplotlib.pyplot as plt
from stackmaster.core import stack
from scipy.signal import sosfiltfilt, butter

dataroot='./data'
dfile=dataroot+"/stackmaster_testdataset.pk"
d=pickle.load(open(dfile,'rb'))

scale=60
data,dt,lag,d_id=[d["data"],d["dt"],d['lag'],d['id']]
tx=np.arange(-lag,lag+0.5*dt,dt)
extent=[-lag,lag,data.shape[0],0]
dn=data.copy()

sos=butter(4,[0.05,0.5],fs=1/dt,btype="bandpass",output='sos')

stack_method="robust"

for i in range(data.shape[0]):
    dn[i,:]=sosfiltfilt(sos,data[i,:]/np.max(np.abs(data[i,:])))

## plot
plt.figure(figsize=(10,5),facecolor="w")
plt.imshow(dn,extent=extent,cmap="seismic",aspect="auto")

dstack=stack(dn,method=stack_method)
plt.plot(tx,scale*dstack+0.5*data.shape[0],'k',lw=2,label=stack_method)
plt.vlines(0,0,data.shape[0],'k')
plt.xlim([-200,200])
plt.ylim([0,data.shape[0]])
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.title(d_id)
plt.xlabel("time (s)",fontsize=14)
plt.ylabel("order",fontsize=14)
plt.legend(fontsize=12)
plt.show()

You should get the following plot:

See https://github.com/xtyangpsp/StackMaster for tutorials and more detailed descriptions.

Any bugs and ideas are welcome. Please file an issue through GitHub https://github.com/xtyangpsp/StackMaster.

• Yang, X., Bryan, J., Okubo, K., Jiang, C., Clements, T., & Denolle, M. A. (2022). Optimal stacking of noise cross-correlation functions. Geophysical Journal International, 232(3), 1600–1618. https://doi.org/10.1093/gji/ggac410