eGf quality control (QC)

Measuring eGf signal-to-noise ratio and pulse width

Import ObsPy module

In [1]:
from obspy import read
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.geodetics.base import gps2dist_azimuth


from obspy.signal.rotate import rotate_ne_rt
from obspy.signal.rotate import rotate_rt_ne
from obspy.signal.rotate import rotate2zne

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

# obspy version =  1.2.2

Import SciPy, NumPy, matplotlib module

In [2]:
import numpy as np
import scipy as sp
import matplotlib as mpl

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

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.dates as mdates 


from scipy import signal
from scipy import ndimage

from scipy.fftpack import fft, ifft
from scipy.linalg import norm

from scipy.stats import pearsonr
from scipy.stats import spearmanr
from scipy.stats import kendalltau

import sys
import os

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

Font size

In [3]:
# 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 moment rate function (MRF)

MRF file in SAC format has 1kHz sampling (interpolated with SAC). I may includ the interpolation part in here or use "raw" MRF file (most of them are 100Hz sampling data)

In [4]:
#st = read("http://ncedc.org/ftp/outgoing/taira/eGf/CLC.stack.sac")

st = read("http://ncedc.org/ftp/outgoing/taira/eGf/stack_indep0.001.sac")
In [5]:
#print(st)

plot MRF

In [6]:
_plot = st.plot(size=(1000,300))
In [7]:
#print(st[0].stats)
In [8]:
tr = st[0]
In [9]:
tr.stats
Out[9]:
         network: 
         station: CLC
        location: 
         channel: 
       starttime: 1970-01-01T00:00:00.000000Z
         endtime: 1970-01-01T00:01:00.008004Z
   sampling_rate: 999.9999389648438
           delta: 0.00100000006103516
            npts: 60009
           calib: 1.0
         _format: SAC
             sac: AttribDict({'delta': 0.001, 'depmin': -2.0186388, 'depmax': 5.8416824, 'b': 0.0, 'e': 60.008003, 't1': 20.042191, 't6': 20.604887, 'stla': 35.815739, 'stlo': -117.59751, 'evla': 35.716, 'evlo': -117.56, 'evdp': 1.92, 'mag': 4.5799999, 'dist': 11.574831, 'az': 342.96835, 'baz': 162.94652, 'gcarc': 0.1040962, 'depmen': 0.0073294407, 'nvhdr': 6, 'npts': 60009, 'iftype': 1, 'leven': 1, 'lpspol': 0, 'lovrok': 1, 'lcalda': 1, 'unused23': 0, 'kstnm': 'CLC', 'kevnm': '38443719'})

get station code

In [10]:
sta = tr.stats.station
print("# sta = ", sta)

# sta =  CLC
In [11]:
#print(tr.data)

get sampling rate

In [12]:
sampling_rateOUT  = int(np.round(tr.stats.sampling_rate))
print("# sampling_rateOUT = ", sampling_rateOUT)

# sampling_rateOUT =  1000
In [13]:
#plt.plot(tr.times(), tr.data)

set parameters

tw is the length of window, st_signal = starting time of signal, st_noise = starting time of noise

In [14]:
delay = 20.0 # sec where we will see MRF onset
tw = 5.0 # sec　window lengh of signal and noise.

# signal part 19.0-24.0
st_signal = 19.0 # sec
et_signal = st_signal + tw # sec

# noise part 5.0-10.sec
st_noise = 5.0 # sec
et_noise = st_noise + tw # sec
In [15]:
time = tr.times() # 
signal_time_array = time[0:int(tw*sampling_rateOUT )] 
noise_time_array = time[0:int(tw*sampling_rateOUT )]
In [16]:
#st_signal / sampling_rateOUT

get signal part in array

In [17]:
st_signal_array = int(sampling_rateOUT * st_signal)
et_signal_array = int(sampling_rateOUT * et_signal)

print("# st_signal_array = ", st_signal_array)
print("# et_signal_array = ", et_signal_array)

# st_signal_array =  19000
# et_signal_array =  24000

get noise part in array

In [18]:
st_noise_array = int(sampling_rateOUT * st_noise)
et_noise_array = int(sampling_rateOUT * et_noise)

print("# st_noise_array = ", st_noise_array)
print("# et_noise_array = ", et_noise_array)

# st_noise_array =  5000
# et_noise_array =  10000

get signal and noise parts

In [19]:
signal = tr.data[st_signal_array:et_signal_array]
noise = tr.data[st_noise_array:et_noise_array]
#print(len(signal))
#print(len(noise))

get signal & noise amplitudes and compute amplitude ratio

In [20]:
signal_amp = norm(signal, 1)/len(signal)
noise_amp = norm(noise, 1)/len(noise)
amp_ratio = signal_amp / noise_amp
print("# amp_ratio = ", amp_ratio)
amp_ratio_out = "{:.3f}".format(amp_ratio)
print("# amp_ratio_out = ", amp_ratio_out)

# amp_ratio =  8.21062243348
# amp_ratio_out =  8.211

get min and max values and their positions in array

In [21]:
max_pos = np.argmax(signal)
print("# max_pos = ", max_pos)
min_pos = np.argmin(signal)
print("# min_pos = ", min_pos)

# max_pos =  1260
# min_pos =  1840
In [22]:
#print(time)

get zero crossing positions in signal array

In [23]:
zero_crossings = np.where(np.diff(np.sign(signal)))[0]
In [24]:
#print(zero_crossings)

get first zero-crossing position after the max value in signal array

This would be the end of MRF pulse

In [25]:
post_cross = zero_crossings[np.where(zero_crossings > max_pos)][0]
print("# post_cross = ", post_cross)

# post_cross =  1605

get first zero-crossing point before the max value in array

This would be the start of MRF pulse

In [26]:
pre_cross = zero_crossings[np.where(zero_crossings <  max_pos)][-1]
print("# pre_cross = ", pre_cross)

# pre_cross =  1042

compute pulse width

In [27]:
pulse_width = (post_cross - pre_cross)/sampling_rateOUT
print("# pulse_width = ", pulse_width, "sec")

pulse_width_out = "{:.3f}".format(pulse_width)
print("# pulse_width_out = ", pulse_width_out, "sec")

# pulse_width =  0.563 sec
# pulse_width_out =  0.563 sec

plot MRF with pulse width

In [28]:
ymin = signal[min_pos] * 1.15
ymax = signal[max_pos] * 1.15

plt.figure(figsize=(18, 4))
plt.plot(signal_time_array, signal, label="signal")
plt.plot(noise_time_array, noise, label="noise")

plt.vlines(x=pre_cross/sampling_rateOUT , ymin=ymin, ymax=ymax, color="black", linewidth = 1.25, linestyle = "--")
plt.vlines(x=post_cross/sampling_rateOUT , ymin=ymin, ymax=ymax, color="black", linewidth = 1.25, linestyle = "--")
plt.xlim([0, tw])
plt.ylim([ymin, ymax])

plt.title(sta+" amp_ratio = "+amp_ratio_out+" pulse_width = "+pulse_width_out+"sec")

plt.ylabel("MRF amplitude")
plt.xlabel("Elapsed time (sec)")
plt.legend(loc = "upper right")
Out[28]:
<matplotlib.legend.Legend at 0x7fac71123520>
In [ ]:
In [ ]:
In [ ]:
In [ ]: