In [573]:
matplotlib inline
In [574]:
from obspy.core import read
from obspy import UTCDateTime
#from obspy.xseed import Parser
#from obspy.signal.rotate import rotate2ZNE
from obspy.core import Stream, Trace, read, AttribDict
from obspy.signal.invsim import paz_2_amplitude_value_of_freq_resp
In [575]:
import matplotlib.pyplot as plt
import pandas as pd
In [576]:
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
#plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 18, 6
In [577]:
import numpy as np
import pandas as pd
from obspy import UTCDateTime
import matplotlib.pyplot as plt
import math
from scipy.optimize import curve_fit
import numpy.random as npr
import os.path
In [578]:
#print(np.__version__)
In [579]:
import obspy as op
#print(op.__version__)
In [580]:
import sys
#print(sys.version)
In [581]:
import warnings
warnings.filterwarnings('ignore')
In [582]:
from IPython.display import HTML
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[582]:
In [583]:
hnm_pdf =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/HNM.dat", sep=" ",header=None)
In [584]:
lnm_pdf =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/LNM.dat", sep=" ",header=None)
In [585]:
data_set_name = "run_071919to091719"
In [586]:
data_set_name =  "2019-07-16T00:00:00.2019-09-18T00:00:00"
In [587]:
#data_set_name =  "2019-10-16T00:00:00.2019-10-23T00:00:00"
In [588]:
BKSOPT = 0
BK63OPT = 0
BK64OPT = 0
BK65OPT = 0
BK69OPT = 0
BK70OPT = 0
BK71OPT = 0
BK72OPT = 0
BK80OPT = 0
BK82OPT = 0
In [589]:
if data_set_name == "2019-07-16T00:00:00.2019-09-18T00:00:00":
    BKSOPT = 1
    BK63OPT = 1
    BK64OPT = 1
    BK65OPT = 1
    BK69OPT = 1
    BK71OPT = 1
    BK72OPT = 1
    BK80OPT = 1
    BK82OPT = 1
In [590]:
if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
    BKSOPT = 1
    BK63OPT = 1
    BK64OPT = 1
    BK65OPT = 1
    BK69OPT = 1
    BK80OPT =1
In [591]:
#data_set_name
In [592]:
#BK72OPT
In [593]:
#if data_set_name != "2019-10-16T00:00:00.2019-10-23T00:00:00":
if data_set_name != "X":
    bb_100s =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/bktest_BB_0.009413_.out", 
                          sep=" ",names=["sta", "zval_med", "nval_med", "eval_med", "nz_val_med", "ez_val_med", "en_val_med"],
                          header=None, skiprows=1)
In [594]:
#bb_200s =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/bktest_BB_0.005132_.out", 
#                          sep=" ",names=["sta", "zval_med", "nval_med", "eval_med", "nz_val_med", "ez_val_med", "en_val_med"],
#                          header=None, skiprows=1)
In [595]:
#if data_set_name != "2019-10-16T00:00:00.2019-10-23T00:00:00":
if data_set_name != "X":
    bb_50s =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/bktest_BB_0.020530_.out", 
                          sep=" ",names=["sta", "zval_med", "nval_med", "eval_med", "nz_val_med", "ez_val_med", "en_val_med"],
                          header=None, skiprows=1)
In [596]:
#if data_set_name != "2019-10-16T00:00:00.2019-10-23T00:00:00":
if data_set_name != "X":
    bb_1s =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/bktest_BB_1.013156_.out", 
                          sep=" ",names=["sta", "zval_med", "nval_med", "eval_med", "nz_val_med", "ez_val_med", "en_val_med"],
                          header=None, skiprows=1)
In [597]:
if BKSOPT == 1:
    bks_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BKS.BK.HHZ.00_.out", 
                          sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                          header=None, skiprows=1)
    bks_hhz_psd_nona = bks_hhz_psd.dropna()
    bks_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BKS.BK.HHN.00_.out", 
                          sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                          header=None, skiprows=1)
    bks_hhn_psd_nona = bks_hhn_psd.dropna()
    bks_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BKS.BK.HHE.00_.out", 
                          sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                          header=None, skiprows=1)
    bks_hhe_psd_nona = bks_hhe_psd.dropna()
In [ ]:
In [598]:
if BK82OPT == 1:
    bk82_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK82.BK.HHZ.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk82_hhz_psd_nona = bk82_hhz_psd.dropna()
    bk82_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK82.BK.HHN.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk82_hhn_psd_nona = bk82_hhn_psd.dropna()
    bk82_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK82.BK.HHE.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                                header=None, skiprows=1)
    bk82_hhe_psd_nona = bk82_hhe_psd.dropna()
In [599]:
if BK80OPT == 1:
    bk80_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK80.BK.HHZ.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk80_hhz_psd_nona = bk80_hhz_psd.dropna()
    bk80_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK80.BK.HHN.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk80_hhn_psd_nona = bk80_hhn_psd.dropna()
    bk80_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK80.BK.HHE.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk80_hhe_psd_nona = bk80_hhe_psd.dropna()
In [600]:
if BK69OPT == 1:
    bk69_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK69.BK.HHZ.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk69_hhz_psd_nona = bk69_hhz_psd.dropna()
    bk69_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK69.BK.HHN.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk69_hhn_psd_nona = bk69_hhn_psd.dropna()
    bk69_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK69.BK.HHE.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk69_hhe_psd_nona = bk69_hhe_psd.dropna()
In [601]:
if BK70OPT == 1:
    bk70_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK70.BK.HHZ.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk70_hhz_psd_nona = bk70_hhz_psd.dropna()
    bk70_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK70.BK.HHN.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk70_hhn_psd_nona = bk70_hhn_psd.dropna()
    bk70_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK70.BK.HHE.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk70_hhe_psd_nona = bk70_hhe_psd.dropna()
In [602]:
if BK71OPT == 1:
    bk71_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK71.BK.HHZ.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk71_hhz_psd_nona = bk71_hhz_psd.dropna()
    bk71_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK71.BK.HHN.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk71_hhn_psd_nona = bk71_hhn_psd.dropna()
    bk71_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK71.BK.HHE.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                                header=None, skiprows=1)
    bk71_hhe_psd_nona = bk71_hhe_psd.dropna()
In [603]:
if BK72OPT == 1:
    bk72_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK72.BK.HHZ.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk72_hhz_psd_nona = bk72_hhz_psd.dropna()
    bk72_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK72.BK.HHN.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk72_hhn_psd_nona = bk72_hhn_psd.dropna()
    bk72_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK72.BK.HHE.00_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk72_hhe_psd_nona = bk72_hhe_psd.dropna()
In [604]:
if BK63OPT == 1:
    bk63_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK63.BK.HHZ.--_.out", 
                                sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk63_hhz_psd_nona = bk63_hhz_psd.dropna()
    bk63_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK63.BK.HHN.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk63_hhn_psd_nona = bk63_hhn_psd.dropna()
    bk63_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK63.BK.HHE.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk63_hhe_psd_nona = bk63_hhe_psd.dropna()
In [605]:
if BK64OPT == 1:
    bk64_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK64.BK.HHZ.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk64_hhz_psd_nona = bk64_hhz_psd.dropna()
    bk64_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK64.BK.HHN.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk64_hhn_psd_nona = bk64_hhn_psd.dropna()
    bk64_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK64.BK.HHE.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk64_hhe_psd_nona = bk64_hhe_psd.dropna()
In [606]:
if BK65OPT == 1:
    bk65_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK65.BK.HHZ.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk65_hhz_psd_nona = bk65_hhz_psd.dropna()
    bk65_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK65.BK.HHN.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk65_hhn_psd_nona = bk65_hhn_psd.dropna()
    bk65_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/"+data_set_name+"/freq.stack_BK65.BK.HHE.--_.out", 
                              sep=" ",names=["freq", "med", "mean", "mode", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
                              header=None, skiprows=1)
    bk65_hhe_psd_nona = bk65_hhe_psd.dropna()
In [ ]:
In [607]:
if BKSOPT == 1:
    ref_hhz_pasd_nona = bks_hhz_psd_nona.copy()
    ref_hhn_pasd_nona = bks_hhn_psd_nona.copy()
    ref_hhe_pasd_nona = bks_hhe_psd_nona.copy()
In [608]:
if BK63OPT == 1:
    bk63_hhz_psd_nona['rmode'] = bk63_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk63_hhn_psd_nona['rmode'] = bk63_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk63_hhe_psd_nona['rmode'] = bk63_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [609]:
if BK64OPT == 1:
    bk64_hhz_psd_nona['rmode'] = bk64_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk64_hhn_psd_nona['rmode'] = bk64_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk64_hhe_psd_nona['rmode'] = bk64_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [610]:
if BK65OPT == 1:
    bk65_hhz_psd_nona['rmode'] = bk65_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk65_hhn_psd_nona['rmode'] = bk65_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk65_hhe_psd_nona['rmode'] = bk65_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [611]:
if BK69OPT == 1:
    bk69_hhz_psd_nona['rmode'] = bk69_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk69_hhn_psd_nona['rmode'] = bk69_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk69_hhe_psd_nona['rmode'] = bk69_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [612]:
if BK70OPT == 1:
    bk70_hhz_psd_nona['rmode'] = bk70_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk70_hhn_psd_nona['rmode'] = bk70_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk70_hhe_psd_nona['rmode'] = bk70_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [613]:
if BK71OPT == 1:
    bk71_hhz_psd_nona['rmode'] = bk71_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk71_hhn_psd_nona['rmode'] = bk71_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk71_hhe_psd_nona['rmode'] = bk71_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [614]:
if BK72OPT == 1:
    bk72_hhz_psd_nona['rmode'] = bk72_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk72_hhn_psd_nona['rmode'] = bk72_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk72_hhe_psd_nona['rmode'] = bk72_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [615]:
if BK80OPT == 1:
    bk80_hhz_psd_nona['rmode'] = bk80_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk80_hhn_psd_nona['rmode'] = bk80_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk80_hhe_psd_nona['rmode'] = bk80_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [616]:
if BK82OPT == 1:
    bk82_hhz_psd_nona['rmode'] = bk82_hhz_psd_nona['mode'] - ref_hhz_pasd_nona['mode']
    bk82_hhn_psd_nona['rmode'] = bk82_hhn_psd_nona['mode'] - ref_hhn_pasd_nona['mode']
    bk82_hhe_psd_nona['rmode'] = bk82_hhe_psd_nona['mode'] - ref_hhe_pasd_nona['mode']
In [617]:
#BKSOPT
In [618]:
print(data_set_name)

2019-07-16T00:00:00.2019-09-18T00:00:00
In [619]:
plt.rcParams['figure.figsize'] = 18, 8
plt.rcParams['figure.figsize'] = 18, 10

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots()
ticks = [10,20,50,100,200,500,1000]
ax.set_xticks(ticks, minor=False)
cmap = plt.get_cmap("tab20c") 

#plt.plot(1/bk66_hhz_psd_nona.freq, bk66_hhz_psd_nona['mode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
#plt.plot(1/bk67_hhz_psd_nona.freq, bk67_hhz_psd_nona['mode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
#plt.plot(1/bk68_hhz_psd_nona.freq, bk68_hhz_psd_nona['mode'], color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')


if BKSOPT == 1:
    plt.plot(1/bks_hhz_psd_nona.freq, bks_hhz_psd_nona['mode'], color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS STS-1')

if BK63OPT == 1:
    plt.plot(1/bk63_hhz_psd_nona.freq, bk63_hhz_psd_nona['mode'], color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK63 Trillium Compact')

if BK64OPT == 1:
    bk64label='BK64 Trillium 120 PH'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk64label='BK64 Trillium 120 PH (glassbeads 3/5mm mix)'
            
    plt.plot(1/bk64_hhz_psd_nona.freq, bk64_hhz_psd_nona['mode'], color=cmap(2), linestyle='solid',linewidth = 3.0, label=bk64label)

if BK65OPT == 1:
    plt.plot(1/bk65_hhz_psd_nona.freq, bk65_hhz_psd_nona['mode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK65 Trillium Horizon 120')

if BK69OPT == 1:
    bk69label='BK69 MBB-2'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk69label='BK69 MBB-2 (sand)'
             
    plt.plot(1/bk69_hhz_psd_nona.freq, bk69_hhz_psd_nona['mode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label=bk69label)

if BK70OPT == 1:
    plt.plot(1/bk70_hhz_psd_nona.freq, bk70_hhz_psd_nona['mode'], color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK70 MBB-2')

if BK71OPT == 1:
    plt.plot(1/bk71_hhz_psd_nona.freq, bk71_hhz_psd_nona['mode'], color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK71 MBB-2')

if BK72OPT == 1:
    plt.plot(1/bk72_hhz_psd_nona.freq, bk72_hhz_psd_nona['mode'], color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK72 MBB-2')

if BK80OPT == 1:
    plt.plot(1/bk80_hhz_psd_nona.freq, bk80_hhz_psd_nona['mode'], color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80 CMG-3T (sand small beads)')

if BK82OPT == 1:
    plt.plot(1/bk82_hhz_psd_nona.freq, bk82_hhz_psd_nona['mode'], color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK82 CMG-3T (sand 3/5mm mix)')

    
plt.plot(hnm_pdf[0], hnm_pdf[1], color='black',linestyle='dashed',linewidth = 2.0,label='HNM')
plt.plot(lnm_pdf[0], lnm_pdf[1], color='black',linestyle='solid',linewidth = 2.0,label='LNM')
plt.xscale("log")
plt.xlim(0.01,10000)
plt.xlim(0.01,100000)
plt.xlim(0.01,2000)

plt.ylim(-190,-80)
plt.ylim(-190,-40)

plt.grid(which='major',color='black',linestyle='-',linewidth = 0.5)
plt.grid(which='minor',color='black',linestyle='-',linewidth = 0.25)

plt.xlabel("Periods (s)", fontsize=16)
plt.ylabel("PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault PSD Com=HHZ "+data_set_name, fontsize=16)

plt.tick_params(labelsize=14)
#plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000,5000,10000], 
#           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000","5000","10000"])
plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000], 
           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000"])

#plt.legend([p3, p2, p1, p4, p5, p6, p7], ["BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)", "USGS-Crest-Mean (HHN,HHE)","BDSN-STS1-Mean (LHN,LHE)"], loc="upper right", fontsize=14)
#plt.legend([p33, p22, p3, p2, p1, p4, p5], ["BK.OAKV.LHN", "BK.OAKV.LHE", "BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)"], loc="upper right", fontsize=14)

plt.legend(loc="upper left", fontsize=14) # 凡例を表示

print("Figure 1: PSD of broadband data in the vertical component")

Figure 1: PSD of broadband data in the vertical component
In [ ]:
In [620]:
plt.rcParams['figure.figsize'] = 18, 10

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots()
ticks = [10,20,50,100,200,500,1000]
ax.set_xticks(ticks, minor=False)
cmap = plt.get_cmap("tab20c") 




if BKSOPT == 1:
    plt.plot(1/bks_hhn_psd_nona.freq, bks_hhn_psd_nona['mode'], color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS STS-1')

if BK63OPT == 1:
    plt.plot(1/bk63_hhn_psd_nona.freq, bk63_hhn_psd_nona['mode'], color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK63 Trillium Compact')

if BK64OPT == 1:
    bk64label='BK64 Trillium 120 PH'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk64label='BK64 Trillium 120 PH (glassbeads 3/5mm mix)'
            
    plt.plot(1/bk64_hhn_psd_nona.freq, bk64_hhn_psd_nona['mode'], color=cmap(2), linestyle='solid',linewidth = 3.0, label=bk64label)

if BK65OPT == 1:
    plt.plot(1/bk65_hhn_psd_nona.freq, bk65_hhn_psd_nona['mode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK65 Trillium Horizon 120')

if BK69OPT == 1:
    bk69label='BK69 MBB-2'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk69label='BK69 MBB-2 (sand)'
             
    plt.plot(1/bk69_hhn_psd_nona.freq, bk69_hhn_psd_nona['mode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label=bk69label)

if BK70OPT == 1:
    plt.plot(1/bk70_hhn_psd_nona.freq, bk70_hhn_psd_nona['mode'], color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK70 MBB-2')

if BK71OPT == 1:
    plt.plot(1/bk71_hhn_psd_nona.freq, bk71_hhn_psd_nona['mode'], color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK71 MBB-2')

if BK72OPT == 1:
    plt.plot(1/bk72_hhn_psd_nona.freq, bk72_hhn_psd_nona['mode'], color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK72 MBB-2')

if BK80OPT == 1:
    plt.plot(1/bk80_hhn_psd_nona.freq, bk80_hhn_psd_nona['mode'], color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80 CMG-3T (sand small beads)')

if BK82OPT == 1:
    plt.plot(1/bk82_hhn_psd_nona.freq, bk82_hhn_psd_nona['mode'], color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK82 CMG-3T (sand 3/5mm mix)')


plt.plot(hnm_pdf[0], hnm_pdf[1], color='black',linestyle='dashed',linewidth = 2.0,label='HNM')
plt.plot(lnm_pdf[0], lnm_pdf[1], color='black',linestyle='solid',linewidth = 2.0,label='LNM')
plt.xscale("log")
plt.xlim(0.01,10000)
plt.xlim(0.01,100000)
plt.xlim(0.01,2000)

plt.ylim(-190,-80)
plt.ylim(-190,-40)

plt.grid(which='major',color='black',linestyle='-',linewidth = 0.5)
plt.grid(which='minor',color='black',linestyle='-',linewidth = 0.25)

plt.xlabel("Periods (s)", fontsize=16)
plt.ylabel("PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault PSD Com=HHN "+data_set_name,fontsize=16)

plt.tick_params(labelsize=14)
#plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000,5000,10000], 
#           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000","5000","10000"])
plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000], 
           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000"])

#plt.legend([p3, p2, p1, p4, p5, p6, p7], ["BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)", "USGS-Crest-Mean (HHN,HHE)","BDSN-STS1-Mean (LHN,LHE)"], loc="upper right", fontsize=14)
#plt.legend([p33, p22, p3, p2, p1, p4, p5], ["BK.OAKV.LHN", "BK.OAKV.LHE", "BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)"], loc="upper right", fontsize=14)

plt.legend(loc="upper left", fontsize=14) # 凡例を表示

print("Figure 2: PSD of broadband data in the N-S component")

Figure 2: PSD of broadband data in the N-S component
In [ ]:
In [621]:
plt.rcParams['figure.figsize'] = 18, 10
#plt.rcParams['figure.figsize'] = 14, 6

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots()
ticks = [10,20,50,100,200,500,1000]
ax.set_xticks(ticks, minor=False)
cmap = plt.get_cmap("tab20c") 



if BKSOPT == 1:
    plt.plot(1/bks_hhe_psd_nona.freq, bks_hhe_psd_nona['mode'], color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS STS-1')

if BK63OPT == 1:
    plt.plot(1/bk63_hhe_psd_nona.freq, bk63_hhe_psd_nona['mode'], color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK63 Trillium Compact')

if BK64OPT == 1:
    bk64label='BK64 Trillium 120 PH'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk64label='BK64 Trillium 120 PH (glassbeads 3/5mm mix)'
            
    plt.plot(1/bk64_hhe_psd_nona.freq, bk64_hhe_psd_nona['mode'], color=cmap(2), linestyle='solid',linewidth = 3.0, label=bk64label)

if BK65OPT == 1:
    plt.plot(1/bk65_hhe_psd_nona.freq, bk65_hhe_psd_nona['mode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK65 Trillium Horizon 120')

if BK69OPT == 1:
    bk69label='BK69 MBB-2'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk69label='BK69 MBB-2 (sand)'
             
    plt.plot(1/bk69_hhe_psd_nona.freq, bk69_hhe_psd_nona['mode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label=bk69label)

if BK70OPT == 1:
    plt.plot(1/bk70_hhe_psd_nona.freq, bk70_hhe_psd_nona['mode'], color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK70 MBB-2')

if BK71OPT == 1:
    plt.plot(1/bk71_hhe_psd_nona.freq, bk71_hhe_psd_nona['mode'], color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK71 MBB-2')

if BK72OPT == 1:
    plt.plot(1/bk72_hhe_psd_nona.freq, bk72_hhe_psd_nona['mode'], color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK72 MBB-2')

if BK80OPT == 1:
    plt.plot(1/bk80_hhe_psd_nona.freq, bk80_hhe_psd_nona['mode'], color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80 CMG-3T (sand small beads)')

if BK82OPT == 1:
    plt.plot(1/bk82_hhe_psd_nona.freq, bk82_hhe_psd_nona['mode'], color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK82 CMG-3T (sand 3/5mm mix)')



plt.plot(hnm_pdf[0], hnm_pdf[1], color='black',linestyle='dashed',linewidth = 2.0,label='HNM')
plt.plot(lnm_pdf[0], lnm_pdf[1], color='black',linestyle='solid',linewidth = 2.0,label='LNM')
plt.xscale("log")

#plt.xlim(0.01,10000)
#plt.xlim(0.01,100000)
plt.xlim(0.01,2000)

plt.ylim(-190,-80)
plt.ylim(-190,-40)

plt.grid(which='major',color='black',linestyle='-',linewidth = 0.5)
plt.grid(which='minor',color='black',linestyle='-',linewidth = 0.25)

plt.xlabel("Periods (s)", fontsize=16)
plt.ylabel("PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault PSD Com=HHE "+data_set_name,fontsize=16)

plt.tick_params(labelsize=14)
#plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000,5000,10000], 
#           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000","5000","10000"])
plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000], 
           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000"])

#plt.legend([p3, p2, p1, p4, p5, p6, p7], ["BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)", "USGS-Crest-Mean (hhe,HHE)","BDSN-STS1-Mean (LHN,LHE)"], loc="upper right", fontsize=14)
#plt.legend([p33, p22, p3, p2, p1, p4, p5], ["BK.OAKV.LHN", "BK.OAKV.LHE", "BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)"], loc="upper right", fontsize=14)

plt.legend(loc="upper left", fontsize=14) # 凡例を表示

print("Figure 3: PSD of broadband data in the E-W component")

Figure 3: PSD of broadband data in the E-W component
In [ ]:
In [622]:
plt.rcParams['figure.figsize'] = 14, 6
plt.rcParams['figure.figsize'] = 18, 10

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots()
ticks = [10,20,50,100,200,500,1000]
ax.set_xticks(ticks, minor=False)
cmap = plt.get_cmap("tab20c") 


#plt.plot(1/bks_hhz_psd_nona.freq, bks_hhz_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')

if BK63OPT == 1:
    plt.plot(1/bk63_hhz_psd_nona.freq, bk63_hhz_psd_nona['rmode'], color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK63 Trillium Compact')

if BK64OPT == 1:
    bk64label='BK64 Trillium 120 PH'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk64label='BK64 Trillium 120 PH (glassbeads 3/5mm mix)'
            
    plt.plot(1/bk64_hhz_psd_nona.freq, bk64_hhz_psd_nona['rmode'], color=cmap(2), linestyle='solid',linewidth = 3.0, label=bk64label)

if BK65OPT == 1:
    plt.plot(1/bk65_hhz_psd_nona.freq, bk65_hhz_psd_nona['rmode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK65 Trillium Horizon 120')

if BK69OPT == 1:
    bk69label='BK69 MBB-2'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk69label='BK69 MBB-2 (sand)'
 
    plt.plot(1/bk69_hhz_psd_nona.freq, bk69_hhz_psd_nona['rmode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label=bk69label)

if BK71OPT == 1:
    plt.plot(1/bk71_hhz_psd_nona.freq, bk71_hhz_psd_nona['rmode'], color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK71 MBB-2')

if BK72OPT == 1:
    plt.plot(1/bk72_hhz_psd_nona.freq, bk72_hhz_psd_nona['rmode'], color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK72 MBB-2')

if BK80OPT == 1:
    plt.plot(1/bk80_hhz_psd_nona.freq, bk80_hhz_psd_nona['rmode'], color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80 CMG-3T (sand small beads)')

if BK82OPT == 1:
    plt.plot(1/bk82_hhz_psd_nona.freq, bk82_hhz_psd_nona['rmode'], color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK82 CMG-3T (sand 3/5mm mix)')


#plt.plot(hnm_pdf[0], hnm_pdf[1], color='black',linestyle='dashed',linewidth = 2.0,label='HNM')
#plt.plot(lnm_pdf[0], lnm_pdf[1], color='black',linestyle='solid',linewidth = 2.0,label='LNM')
plt.xscale("log")
plt.xlim(0.01,10000)
plt.xlim(0.01,2000)
plt.ylim(-20,80)

plt.grid(which='major',color='black',linestyle='-',linewidth = 0.5)
plt.grid(which='minor',color='black',linestyle='-',linewidth = 0.25)

plt.xlabel("Periods (s)", fontsize=16)
plt.ylabel("Relative PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault rPSD (ref:BKS) Com=HHZ "+data_set_name,fontsize=16)

plt.tick_params(labelsize=14)
#plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000,5000,10000], 
#           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000","5000","10000"])
plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000], 
           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000"])

#plt.legend([p3, p2, p1, p4, p5, p6, p7], ["BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)", "USGS-Crest-Mean (HHN,HHE)","BDSN-STS1-Mean (LHN,LHE)"], loc="upper right", fontsize=14)
#plt.legend([p33, p22, p3, p2, p1, p4, p5], ["BK.OAKV.LHN", "BK.OAKV.LHE", "BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)"], loc="upper right", fontsize=14)

plt.legend(loc="upper left", fontsize=14) # 凡例を表示

print("Figure 4: Relative PSD (BKS as reference) of broadband data in the vertical component")

Figure 4: Relative PSD (BKS as reference) of broadband data in the vertical component
In [ ]:
In [ ]:
In [623]:
plt.rcParams['figure.figsize'] = 14, 6
plt.rcParams['figure.figsize'] = 18, 10

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots()
ticks = [10,20,50,100,200,500,1000]
ax.set_xticks(ticks, minor=False)
cmap = plt.get_cmap("tab20c") 


#plt.plot(1/bks_hhn_psd_nona.freq, bks_hhn_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')


if BK63OPT == 1:
    plt.plot(1/bk63_hhn_psd_nona.freq, bk63_hhn_psd_nona['rmode'], color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK63 Trillium Compact')

if BK64OPT == 1:
    bk64label='BK64 Trillium 120 PH'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk64label='BK64 Trillium 120 PH (glassbeads 3/5mm mix)'
            
    plt.plot(1/bk64_hhn_psd_nona.freq, bk64_hhn_psd_nona['rmode'], color=cmap(2), linestyle='solid',linewidth = 3.0, label=bk64label)

if BK65OPT == 1:
    plt.plot(1/bk65_hhn_psd_nona.freq, bk65_hhn_psd_nona['rmode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK65 Trillium Horizon 120')

if BK69OPT == 1:
    bk69label='BK69 MBB-2'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk69label='BK69 MBB-2 (sand)'
 
    plt.plot(1/bk69_hhn_psd_nona.freq, bk69_hhn_psd_nona['rmode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label=bk69label)

if BK71OPT == 1:
    plt.plot(1/bk71_hhn_psd_nona.freq, bk71_hhn_psd_nona['rmode'], color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK71 MBB-2')

if BK72OPT == 1:
    plt.plot(1/bk72_hhn_psd_nona.freq, bk72_hhn_psd_nona['rmode'], color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK72 MBB-2')

if BK80OPT == 1:
    plt.plot(1/bk80_hhn_psd_nona.freq, bk80_hhn_psd_nona['rmode'], color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80 CMG-3T (sand small beads)')

if BK82OPT == 1:
    plt.plot(1/bk82_hhn_psd_nona.freq, bk82_hhn_psd_nona['rmode'], color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK82 CMG-3T (sand 3/5mm mix)')


#plt.plot(hnm_pdf[0], hnm_pdf[1], color='black',linestyle='dashed',linewidth = 2.0,label='HNM')
#plt.plot(lnm_pdf[0], lnm_pdf[1], color='black',linestyle='solid',linewidth = 2.0,label='LNM')
plt.xscale("log")
plt.xlim(0.01,10000)
plt.xlim(0.01,2000)
plt.ylim(-10,50)
plt.ylim(-20,80)

plt.grid(which='major',color='black',linestyle='-',linewidth = 0.5)
plt.grid(which='minor',color='black',linestyle='-',linewidth = 0.25)

plt.xlabel("Periods (s)", fontsize=16)
plt.ylabel("Relative PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault rPSD (ref:BKS) Com=HHN "+data_set_name,fontsize=16)

plt.tick_params(labelsize=14)
#plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000,5000,10000], 
#           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000","5000","10000"])
plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000], 
           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000"])

#plt.legend([p3, p2, p1, p4, p5, p6, p7], ["BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)", "USGS-Crest-Mean (HHN,HHE)","BDSN-STS1-Mean (LHN,LHE)"], loc="upper right", fontsize=14)
#plt.legend([p33, p22, p3, p2, p1, p4, p5], ["BK.OAKV.LHN", "BK.OAKV.LHE", "BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)"], loc="upper right", fontsize=14)

plt.legend(loc="upper left", fontsize=14) # 凡例を表示

print("Figure 5: Relative PSD (BKS as reference) of broadband data in the N-S component")

Figure 5: Relative PSD (BKS as reference) of broadband data in the N-S component
In [ ]:
In [624]:
plt.rcParams['figure.figsize'] = 14, 6
plt.rcParams['figure.figsize'] = 18, 10

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, ax = plt.subplots()
ticks = [10,20,50,100,200,500,1000]
ax.set_xticks(ticks, minor=False)
cmap = plt.get_cmap("tab20c") 


#plt.plot(1/bks_hhe_psd_nona.freq, bks_hhe_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')



if BK63OPT == 1:
    plt.plot(1/bk63_hhe_psd_nona.freq, bk63_hhe_psd_nona['rmode'], color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK63 Trillium Compact')

if BK64OPT == 1:
    bk64label='BK64 Trillium 120 PH'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk64label='BK64 Trillium 120 PH (glassbeads 3/5mm mix)'
            
    plt.plot(1/bk64_hhe_psd_nona.freq, bk64_hhe_psd_nona['rmode'], color=cmap(2), linestyle='solid',linewidth = 3.0, label=bk64label)

if BK65OPT == 1:
    plt.plot(1/bk65_hhe_psd_nona.freq, bk65_hhe_psd_nona['rmode'], color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK65 Trillium Horizon 120')

if BK69OPT == 1:
    bk69label='BK69 MBB-2'
    if data_set_name == "2019-10-16T00:00:00.2019-10-23T00:00:00":
          bk69label='BK69 MBB-2 (sand)'
 
    plt.plot(1/bk69_hhe_psd_nona.freq, bk69_hhe_psd_nona['rmode'], color=cmap(4), linestyle='solid',linewidth = 3.0, label=bk69label)

if BK71OPT == 1:
    plt.plot(1/bk71_hhe_psd_nona.freq, bk71_hhe_psd_nona['rmode'], color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK71 MBB-2')

if BK72OPT == 1:
    plt.plot(1/bk72_hhe_psd_nona.freq, bk72_hhe_psd_nona['rmode'], color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK72 MBB-2')

if BK80OPT == 1:
    plt.plot(1/bk80_hhe_psd_nona.freq, bk80_hhe_psd_nona['rmode'], color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80 CMG-3T (sand small beads)')

if BK82OPT == 1:
    plt.plot(1/bk82_hhe_psd_nona.freq, bk82_hhe_psd_nona['rmode'], color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK82 CMG-3T (sand 3/5mm mix)')



#plt.plot(hnm_pdf[0], hnm_pdf[1], color='black',linestyle='dashed',linewidth = 2.0,label='HNM')
#plt.plot(lnm_pdf[0], lnm_pdf[1], color='black',linestyle='solid',linewidth = 2.0,label='LNM')
plt.xscale("log")
plt.xlim(0.01,10000)
plt.xlim(0.01,2000)
plt.ylim(-10,50)
plt.ylim(-20,80)
plt.grid(which='major',color='black',linestyle='-',linewidth = 0.5)
plt.grid(which='minor',color='black',linestyle='-',linewidth = 0.25)

plt.xlabel("Periods (s)", fontsize=16)
plt.ylabel("Relative PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault rPSD (ref:BKS) Com=HHE "+data_set_name,fontsize=16)

plt.tick_params(labelsize=14)
#plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000,5000,10000], 
#           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000","5000","10000"])
plt.xticks([0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20, 50,100,200,500,1000,2000], 
           ["0.01","0.02","0.05","0.1","0.2","0.5","1","2","5","10","20","50","100","200","500","1000","2000"])

#plt.legend([p3, p2, p1, p4, p5, p6, p7], ["BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)", "USGS-Crest-Mean (hhe,HHE)","BDSN-STS1-Mean (LHN,LHE)"], loc="upper right", fontsize=14)
#plt.legend([p33, p22, p3, p2, p1, p4, p5], ["BK.OAKV.LHN", "BK.OAKV.LHE", "BK.JEPS.LHN", "BK.JEPS.LHE", "High Noise Model", "TA-Mean (LHN,LHE)", "BDSN-RefTek-Mean (LHN,LHE)"], loc="upper right", fontsize=14)

plt.legend(loc="upper left", fontsize=14) # 凡例を表示

print("Figure 6: Relative PSD (BKS as reference) of broadband data in the E-W component")

Figure 6: Relative PSD (BKS as reference) of broadband data in the E-W component
In [ ]:
In [625]:
#plt.rcParams['figure.figsize'] = 11, 6
plt.rcParams['figure.figsize'] = 15, 5
plt.xlabel('Station',fontsize=16)

sta_list = bb_100s.sta.as_matrix()
zval_med = bb_100s.zval_med.as_matrix()
nval_med = bb_100s.nval_med.as_matrix()
eval_med = bb_100s.eval_med.as_matrix()

x_position = np.arange(len(sta_list))
plt.bar(x_position - 0.3, zval_med, width=0.3, label='HHZ')
plt.bar(x_position, nval_med, width=0.3, label='HHN')                    
plt.bar(x_position + 0.3, eval_med, width=0.3, label='HHE')

plt.legend(fontsize=16)
#plt.xticks(x_position + 0.2, x)
plt.xticks(x_position + 0.1, sta_list)

plt.ylabel('PSD at 100-s period', fontsize=16)
#plt.ylabel('PSD at 200-s period', fontsize=16)
plt.title("Byerly Vault PSDs at 100 s period "+data_set_name,fontsize=16)

# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(-120, -180)
plt.ylim(-120, -190)
plt.legend(loc="upper right", fontsize=14) # 凡例を表示
plt.tick_params(labelsize=14)


print("Figure 7: PSD of broadband data at 100-s period")

plt.show()

Figure 7: PSD of broadband data at 100-s period
In [ ]:
In [ ]:
In [626]:
#plt.rcParams['figure.figsize'] = 11, 6
plt.rcParams['figure.figsize'] = 15, 5
plt.xlabel('Station',fontsize=16)

sta_list = bb_50s.sta.as_matrix()
zval_med = bb_50s.zval_med.as_matrix()
nval_med = bb_50s.nval_med.as_matrix()
eval_med = bb_50s.eval_med.as_matrix()

x_position = np.arange(len(sta_list))
plt.bar(x_position - 0.3, zval_med, width=0.3, label='HHZ')
plt.bar(x_position, nval_med, width=0.3, label='HHN')                    
plt.bar(x_position + 0.3, eval_med, width=0.3, label='HHE')

plt.legend(fontsize=16)
#plt.xticks(x_position + 0.2, x)
plt.xticks(x_position + 0.1, sta_list)

plt.ylabel('PSD at 50-s period', fontsize=16)
#plt.ylabel('PSD at 200-s period', fontsize=16)
plt.title("Byerly Vault PSDs at 50 s period "+data_set_name,fontsize=16)

# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(-120, -190)
plt.legend(loc="upper right", fontsize=14) # 凡例を表示
plt.tick_params(labelsize=14)


print("Figure 8: PSD of broadband data at 50-s period")


plt.show()

Figure 8: PSD of broadband data at 50-s period
In [627]:
#plt.rcParams['figure.figsize'] = 11, 6
plt.rcParams['figure.figsize'] = 15, 5
plt.xlabel('Station',fontsize=16)

sta_list = bb_1s.sta.as_matrix()
zval_med = bb_1s.zval_med.as_matrix()
nval_med = bb_1s.nval_med.as_matrix()
eval_med = bb_1s.eval_med.as_matrix()

x_position = np.arange(len(sta_list))
plt.bar(x_position - 0.3, zval_med, width=0.3, label='HHZ')
plt.bar(x_position, nval_med, width=0.3, label='HHN')                    
plt.bar(x_position + 0.3, eval_med, width=0.3, label='HHE')

plt.legend(fontsize=16)
#plt.xticks(x_position + 0.2, x)
plt.xticks(x_position + 0.1, sta_list)

plt.ylabel('PSD at 1-s period', fontsize=16)
#plt.ylabel('PSD at 200-s period', fontsize=16)
plt.title("Byerly Vault PSDs at 1 s period "+data_set_name,fontsize=16)

# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(-120, -140)
plt.legend(loc="upper right", fontsize=14) # 凡例を表示
plt.tick_params(labelsize=14)


print("Figure 9: PSD of broadband data at 1-s period")


plt.show()

Figure 9: PSD of broadband data at 1-s period
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: