Stretching waveform

making synthetic waveform by "stretching".

Import ObsPy module

from obspy import read
import obspy as ob

print("# obspy version = ",ob.__version__)

# obspy version =  1.2.2

Import SciPy, NumPy, matplotlib module

import numpy as np
import scipy as sp
import matplotlib as mpl

import matplotlib.pyplot as plt

from scipy import interpolate

print("# numpy version = ",np.__version__)
print("# scipy version = ",sp.__version__)
print("# matplotlib version = ",mpl.__version__)

# numpy version =  1.19.1
# scipy version =  1.5.2
# matplotlib version =  3.3.1

Font size

# font size
SMALL_SIZE = 16
MEDIUM_SIZE = 18
BIGGER_SIZE = 20

#SMALL_SIZE = 32
#MEDIUM_SIZE = 32
#BIGGER_SIZE = 36

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title

Read waveform

this is the original waveform.

st = read("http://ncedc.org/ftp/outgoing/taira/LabEQ/PZT08.AE.EP1.00.D.1970.001.052401.2248.sac")
_plot = st.plot(size=(1000,300))

get "time"

time is np.array

tr = st[0]
time = tr.times() 
time_shift = time.copy() # making "time_shift". this would be modified below. 
# or may be we can simply make time_shift = np.zeros(len(time))
#time_shift = np.zeros(len(time))

#type(time)

#plt.figure(figsize=(18, 4))
#plt.plot(time, tr.data)

stretch_coeff

stretch_coeff = 0.01 
#stretch_coeff = 0.05
#stretch_coeff = 0.1
#stretch_coeff = -0.1

# 0.01 means 0.01sec "delay" at elapse time 1sec -> new "time" is 1.01sec
# positive value -> velocity decrese case. negative value -> velocity increse case.

compute time shift

#debug = 1 # if we want to see all outputs
debug = 0

for i, data in enumerate(tr.data):
    #print("# i = ", i)
    shiftT=i*tr.stats.delta*stretch_coeff
    time_shift[i] = time[i] + shiftT 
    if debug:
        print("# shiftT = ", shiftT, " time_shift = ",time_shift[i], "time = ",time[i] )

plot waveforms

plt.figure(figsize=(18, 4))

plt.plot(time, tr.data, label="original") # original
plt.plot(time_shift, tr.data, label="new") # data with updated time stamp
plt.legend()


plt.xlabel("Time (s)")
plt.ylabel("Amplitude")

Text(0, 0.5, 'Amplitude')

plot lag time between original and new data

This lag time should be ideally obtained from moving window cross-correlation analysis

plt.figure(figsize=(18, 4))

lag_time=time_shift-time
plt.plot(time, time_shift-time)

plt.xlabel("Time (s)")
plt.ylabel("Lag time (s)")

Text(0, 0.5, 'Lag time (s)')

trim data

st_data_point = 0 # start data point
et_data_point = 450 # end data point. as test data has 0.1-sec sampling rate, this would be 45 seconds 

st_time = st_data_point * tr.stats.delta
et_time = et_data_point * tr.stats.delta

print("# st_time (s) = ", st_time, " et_time (s) = ", et_time)

# st_time (s) =  0.0  et_time (s) =  45.0

# trim original data
time_org = time[st_data_point:et_data_point]
amp_org = tr.data[st_data_point:et_data_point]

interpolate "new data"

so that we can have "same" time stamp for new data. we use "interpolate" from scipy

f = interpolate.interp1d(time_shift, tr.data, kind="cubic")
time_new = np.arange(st_time, et_time, tr.stats.delta)
amp_new = f(time_new)

plot waveforms

original and interpolated new data. both should have the same sampling time and the same data length

plt.figure(figsize=(18, 4))
plt.plot(time_org, amp_org, label="original")
plt.plot(time_new, amp_new, label="new")
plt.legend(loc = "upper right")

plt.xlabel("Time (s)")
plt.ylabel("Amplitude")

Text(0, 0.5, 'Amplitude')

save waveform data

original data -> PZT08.mseed and new data -> PZT09.mseed

st1 = st.copy()
st1[0].data = amp_org
st1[0].time = time_org
st1.write("PZT08.mseed", format="MSEED")

st2 = st.copy()

st2[0].data = amp_new
st2[0].time = time_new
st2[0].stats.station = "PZT09"
#print(st2)

st2.write("PZT09.mseed", format="MSEED")