In [312]:
matplotlib inline
In [313]:
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 [314]:
import matplotlib.pyplot as plt
import pandas as pd
In [315]:
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
#plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 11, 6
In [316]:
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 [317]:
#print(np.__version__)
In [318]:
import obspy as op
#print(op.__version__)
In [319]:
import sys
#print(sys.version)
In [320]:
import warnings
warnings.filterwarnings('ignore')
In [321]:
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[321]:
In [322]:
hnm_pdf =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/HNM.dat", sep=" ",header=None)
In [323]:
lnm_pdf =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/LNM.dat", sep=" ",header=None)
In [324]:
bb_100s =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/bktest_BB_0.010265_.out",
sep=" ",names=["sta", "zval_med", "nval_med", "eval_med", "nz_val_med", "ez_val_med", "en_val_med"],
header=None, skiprows=1)
In [325]:
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 [326]:
bb_50s =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/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 [327]:
bk70_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK70.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK70.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK70.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk70_hhe_psd_nona = bk70_hhe_psd.dropna()
In [328]:
bks_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BKS.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BKS.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BKS.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bks_hhe_psd_nona = bks_hhe_psd.dropna()
In [329]:
bk63_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK63.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK63.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK63.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk63_hhe_psd_nona = bk63_hhe_psd.dropna()
In [330]:
bk64_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK64.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK64.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK64.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk64_hhe_psd_nona = bk64_hhe_psd.dropna()
In [331]:
bk65_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK65.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK65.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK65.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk65_hhe_psd_nona = bk65_hhe_psd.dropna()
In [332]:
bk66_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK66.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk66_hhz_psd_nona = bk66_hhz_psd.dropna()
bk66_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK66.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk66_hhn_psd_nona = bk66_hhn_psd.dropna()
bk66_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK66.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk66_hhe_psd_nona = bk66_hhe_psd.dropna()
In [333]:
bk67_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK67.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk67_hhz_psd_nona = bk67_hhz_psd.dropna()
bk67_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK67.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk67_hhn_psd_nona = bk67_hhn_psd.dropna()
bk67_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK67.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk67_hhe_psd_nona = bk67_hhe_psd.dropna()
In [334]:
bk68_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK68.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk68_hhz_psd_nona = bk68_hhz_psd.dropna()
bk68_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK68.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk68_hhn_psd_nona = bk68_hhn_psd.dropna()
bk68_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK68.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk68_hhe_psd_nona = bk68_hhe_psd.dropna()
In [335]:
bk69_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK69.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK69.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK69.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk69_hhe_psd_nona = bk69_hhe_psd.dropna()
In [336]:
bk80_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK80.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK80.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK80.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk80_hhe_psd_nona = bk80_hhe_psd.dropna()
In [337]:
bk81_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK81.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk81_hhz_psd_nona = bk81_hhz_psd.dropna()
bk81_hhn_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK81.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk81_hhn_psd_nona = bk81_hhn_psd.dropna()
bk81_hhe_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK81.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk81_hhe_psd_nona = bk81_hhe_psd.dropna()
In [338]:
bk82_hhz_psd =pd.read_csv("http://ncedc.org/ftp/outgoing/taira/WQC/ByerlyVault_PSD/freq.stack_BK82.BK.HHZ.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK82.BK.HHN.00_.out",
sep=" ",names=["freq", "med", "mean", "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/freq.stack_BK82.BK.HHE.00_.out",
sep=" ",names=["freq", "med", "mean", "l95", "h95", "min", "max", "sd","var","sta","net","com","loc"],
header=None, skiprows=1)
bk82_hhe_psd_nona = bk82_hhe_psd.dropna()
In [339]:
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 [340]:
bk63_hhz_psd_nona.rmed = bk63_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk63_hhn_psd_nona.rmed = bk63_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk63_hhe_psd_nona.rmed = bk63_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [341]:
bk64_hhz_psd_nona.rmed = bk64_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk64_hhn_psd_nona.rmed = bk64_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk64_hhe_psd_nona.rmed = bk64_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [342]:
bk65_hhz_psd_nona.rmed = bk65_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk65_hhn_psd_nona.rmed = bk65_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk65_hhe_psd_nona.rmed = bk65_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [343]:
bk66_hhz_psd_nona.rmed = bk66_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk66_hhn_psd_nona.rmed = bk66_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk66_hhe_psd_nona.rmed = bk66_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [344]:
bk67_hhz_psd_nona.rmed = bk67_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk67_hhn_psd_nona.rmed = bk67_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk67_hhe_psd_nona.rmed = bk67_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [345]:
bk68_hhz_psd_nona.rmed = bk68_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk68_hhn_psd_nona.rmed = bk68_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk68_hhe_psd_nona.rmed = bk68_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [346]:
bk69_hhz_psd_nona.rmed = bk69_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk69_hhn_psd_nona.rmed = bk69_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk69_hhe_psd_nona.rmed = bk69_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [347]:
bk70_hhz_psd_nona.rmed = bk70_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk70_hhn_psd_nona.rmed = bk70_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk70_hhe_psd_nona.rmed = bk70_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [348]:
bk80_hhz_psd_nona.rmed = bk80_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk80_hhn_psd_nona.rmed = bk80_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk80_hhe_psd_nona.rmed = bk80_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [349]:
bk81_hhz_psd_nona.rmed = bk81_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk81_hhn_psd_nona.rmed = bk81_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk81_hhe_psd_nona.rmed = bk81_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [350]:
bk82_hhz_psd_nona.rmed = bk82_hhz_psd_nona.med - ref_hhz_pasd_nona.med
bk82_hhn_psd_nona.rmed = bk82_hhn_psd_nona.med - ref_hhn_pasd_nona.med
bk82_hhe_psd_nona.rmed = bk82_hhe_psd_nona.med - ref_hhe_pasd_nona.med
In [351]:
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")
plt.plot(1/bks_hhz_psd_nona.freq, bks_hhz_psd_nona.med, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')
plt.plot(1/bk63_hhz_psd_nona.freq, bk63_hhz_psd_nona.med, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BK63')
plt.plot(1/bk64_hhz_psd_nona.freq, bk64_hhz_psd_nona.med, color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK64')
plt.plot(1/bk65_hhz_psd_nona.freq, bk65_hhz_psd_nona.med, color=cmap(2), linestyle='solid',linewidth = 3.0, label='BK65')
plt.plot(1/bk66_hhz_psd_nona.freq, bk66_hhz_psd_nona.med, color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
plt.plot(1/bk67_hhz_psd_nona.freq, bk67_hhz_psd_nona.med, color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
plt.plot(1/bk68_hhz_psd_nona.freq, bk68_hhz_psd_nona.med, color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')
plt.plot(1/bk69_hhz_psd_nona.freq, bk69_hhz_psd_nona.med, color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK69')
plt.plot(1/bk70_hhz_psd_nona.freq, bk70_hhz_psd_nona.med, color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK70')
plt.plot(1/bk80_hhz_psd_nona.freq, bk80_hhz_psd_nona.med, color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80')
plt.plot(1/bk81_hhz_psd_nona.freq, bk81_hhz_psd_nona.med, color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK81')
plt.plot(1/bk82_hhz_psd_nona.freq, bk82_hhz_psd_nona.med, color=cmap(10), linestyle='solid',linewidth = 3.0, label='BK82')
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.ylim(-190,-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("PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault PSD Com=HHZ",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.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 right", fontsize=14) # 凡例を表示
print("Figure 1: PSD of broadband data in the vertical component")
In [352]:
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")
plt.plot(1/bks_hhn_psd_nona.freq, bks_hhn_psd_nona.med, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')
plt.plot(1/bk63_hhn_psd_nona.freq, bk63_hhn_psd_nona.med, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BK63')
plt.plot(1/bk64_hhn_psd_nona.freq, bk64_hhn_psd_nona.med, color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK64')
plt.plot(1/bk65_hhn_psd_nona.freq, bk65_hhn_psd_nona.med, color=cmap(2), linestyle='solid',linewidth = 3.0, label='BK65')
plt.plot(1/bk66_hhn_psd_nona.freq, bk66_hhn_psd_nona.med, color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
plt.plot(1/bk67_hhn_psd_nona.freq, bk67_hhn_psd_nona.med, color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
plt.plot(1/bk68_hhn_psd_nona.freq, bk68_hhn_psd_nona.med, color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')
plt.plot(1/bk69_hhn_psd_nona.freq, bk69_hhn_psd_nona.med, color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK69')
plt.plot(1/bk70_hhn_psd_nona.freq, bk70_hhn_psd_nona.med, color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK70')
plt.plot(1/bk80_hhn_psd_nona.freq, bk80_hhn_psd_nona.med, color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80')
plt.plot(1/bk81_hhn_psd_nona.freq, bk81_hhn_psd_nona.med, color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK81')
plt.plot(1/bk82_hhn_psd_nona.freq, bk82_hhn_psd_nona.med, color=cmap(10), linestyle='solid',linewidth = 3.0, label='BK82')
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.ylim(-190,-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("PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault PSD Com=HHN",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.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 right", fontsize=14) # 凡例を表示
print("Figure 2: PSD of broadband data in the N-S component")
In [353]:
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")
plt.plot(1/bks_hhe_psd_nona.freq, bks_hhe_psd_nona.med, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')
plt.plot(1/bk63_hhe_psd_nona.freq, bk63_hhe_psd_nona.med, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BK63')
plt.plot(1/bk64_hhe_psd_nona.freq, bk64_hhe_psd_nona.med, color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK64')
plt.plot(1/bk65_hhe_psd_nona.freq, bk65_hhe_psd_nona.med, color=cmap(2), linestyle='solid',linewidth = 3.0, label='BK65')
plt.plot(1/bk66_hhe_psd_nona.freq, bk66_hhe_psd_nona.med, color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
plt.plot(1/bk67_hhe_psd_nona.freq, bk67_hhe_psd_nona.med, color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
plt.plot(1/bk68_hhe_psd_nona.freq, bk68_hhe_psd_nona.med, color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')
plt.plot(1/bk69_hhe_psd_nona.freq, bk69_hhe_psd_nona.med, color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK69')
plt.plot(1/bk70_hhe_psd_nona.freq, bk70_hhe_psd_nona.med, color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK70')
plt.plot(1/bk80_hhe_psd_nona.freq, bk80_hhe_psd_nona.med, color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80')
plt.plot(1/bk81_hhe_psd_nona.freq, bk81_hhe_psd_nona.med, color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK81')
plt.plot(1/bk82_hhe_psd_nona.freq, bk82_hhe_psd_nona.med, color=cmap(10), linestyle='solid',linewidth = 3.0, label='BK82')
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.ylim(-190,-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("PSD Power [10log(m**2/sec**4/Hz)] (dB)", fontsize=16)
plt.title("Byerly Vault PSD Com=HHE",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.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 right", fontsize=14) # 凡例を表示
print("Figure 3: PSD of broadband data in the E-W component")
In [354]:
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")
#plt.plot(1/bks_hhz_psd_nona.freq, bks_hhz_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')
plt.plot(1/bk63_hhz_psd_nona.freq, bk63_hhz_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BK63')
plt.plot(1/bk64_hhz_psd_nona.freq, bk64_hhz_psd_nona.rmed, color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK64')
plt.plot(1/bk65_hhz_psd_nona.freq, bk65_hhz_psd_nona.rmed, color=cmap(2), linestyle='solid',linewidth = 3.0, label='BK65')
plt.plot(1/bk66_hhz_psd_nona.freq, bk66_hhz_psd_nona.rmed, color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
plt.plot(1/bk67_hhz_psd_nona.freq, bk67_hhz_psd_nona.rmed, color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
plt.plot(1/bk68_hhz_psd_nona.freq, bk68_hhz_psd_nona.rmed, color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')
plt.plot(1/bk69_hhz_psd_nona.freq, bk69_hhz_psd_nona.rmed, color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK69')
plt.plot(1/bk70_hhz_psd_nona.freq, bk70_hhz_psd_nona.rmed, color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK70')
plt.plot(1/bk80_hhz_psd_nona.freq, bk80_hhz_psd_nona.rmed, color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80')
plt.plot(1/bk81_hhz_psd_nona.freq, bk81_hhz_psd_nona.rmed, color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK81')
plt.plot(1/bk82_hhz_psd_nona.freq, bk82_hhz_psd_nona.rmed, color=cmap(10), linestyle='solid',linewidth = 3.0, label='BK82')
#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.ylim(-10,50)
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",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.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 right", fontsize=14) # 凡例を表示
print("Figure 4: Relative PSD (BKS as reference) of broadband data in the vertical component")
In [355]:
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")
#plt.plot(1/bks_hhn_psd_nona.freq, bks_hhn_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')
plt.plot(1/bk63_hhn_psd_nona.freq, bk63_hhn_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BK63')
plt.plot(1/bk64_hhn_psd_nona.freq, bk64_hhn_psd_nona.rmed, color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK64')
plt.plot(1/bk65_hhn_psd_nona.freq, bk65_hhn_psd_nona.rmed, color=cmap(2), linestyle='solid',linewidth = 3.0, label='BK65')
plt.plot(1/bk66_hhn_psd_nona.freq, bk66_hhn_psd_nona.rmed, color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
plt.plot(1/bk67_hhn_psd_nona.freq, bk67_hhn_psd_nona.rmed, color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
plt.plot(1/bk68_hhn_psd_nona.freq, bk68_hhn_psd_nona.rmed, color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')
plt.plot(1/bk69_hhn_psd_nona.freq, bk69_hhn_psd_nona.rmed, color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK69')
plt.plot(1/bk70_hhn_psd_nona.freq, bk70_hhn_psd_nona.rmed, color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK70')
plt.plot(1/bk80_hhn_psd_nona.freq, bk80_hhn_psd_nona.rmed, color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80')
plt.plot(1/bk81_hhn_psd_nona.freq, bk81_hhn_psd_nona.rmed, color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK81')
plt.plot(1/bk82_hhn_psd_nona.freq, bk82_hhn_psd_nona.rmed, color=cmap(10), linestyle='solid',linewidth = 3.0, label='BK82')
#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.ylim(-10,50)
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",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.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 right", fontsize=14) # 凡例を表示
print("Figure 5: Relative PSD (BKS as reference) of broadband data in the N-S component")
In [356]:
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")
#plt.plot(1/bks_hhe_psd_nona.freq, bks_hhe_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BKS')
plt.plot(1/bk63_hhe_psd_nona.freq, bk63_hhe_psd_nona.rmed, color=cmap(0), linestyle='solid',linewidth = 3.0, label='BK63')
plt.plot(1/bk64_hhe_psd_nona.freq, bk64_hhe_psd_nona.rmed, color=cmap(1), linestyle='solid',linewidth = 3.0, label='BK64')
plt.plot(1/bk65_hhe_psd_nona.freq, bk65_hhe_psd_nona.rmed, color=cmap(2), linestyle='solid',linewidth = 3.0, label='BK65')
plt.plot(1/bk66_hhe_psd_nona.freq, bk66_hhe_psd_nona.rmed, color=cmap(3), linestyle='solid',linewidth = 3.0, label='BK66')
plt.plot(1/bk67_hhe_psd_nona.freq, bk67_hhe_psd_nona.rmed, color=cmap(4), linestyle='solid',linewidth = 3.0, label='BK67')
plt.plot(1/bk68_hhe_psd_nona.freq, bk68_hhe_psd_nona.rmed, color=cmap(5), linestyle='solid',linewidth = 3.0, label='BK68')
plt.plot(1/bk69_hhe_psd_nona.freq, bk69_hhe_psd_nona.rmed, color=cmap(6), linestyle='solid',linewidth = 3.0, label='BK69')
plt.plot(1/bk70_hhe_psd_nona.freq, bk70_hhe_psd_nona.rmed, color=cmap(7), linestyle='solid',linewidth = 3.0, label='BK70')
plt.plot(1/bk80_hhe_psd_nona.freq, bk80_hhe_psd_nona.rmed, color=cmap(8), linestyle='solid',linewidth = 3.0, label='BK80')
plt.plot(1/bk81_hhe_psd_nona.freq, bk81_hhe_psd_nona.rmed, color=cmap(9), linestyle='solid',linewidth = 3.0, label='BK81')
plt.plot(1/bk82_hhe_psd_nona.freq, bk82_hhe_psd_nona.rmed, color=cmap(10), linestyle='solid',linewidth = 3.0, label='BK82')
#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.ylim(-10,50)
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",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.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 right", fontsize=14) # 凡例を表示
print("Figure 6: Relative PSD (BKS as reference) of broadband data in the E-W component")
In [357]:
#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",fontsize=16)
# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(-120, -180)
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()
In [358]:
#plt.rcParams['figure.figsize'] = 11, 6
plt.rcParams['figure.figsize'] = 15, 5
plt.xlabel('Station',fontsize=16)
sta_list = bb_200s.sta.as_matrix()
zval_med = bb_200s.zval_med.as_matrix()
nval_med = bb_200s.nval_med.as_matrix()
eval_med = bb_200s.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 200-s period', fontsize=16)
#plt.ylabel('PSD at 200-s period', fontsize=16)
plt.title("Byerly Vault PSDs at 200 s period",fontsize=16)
# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(-120, -180)
plt.legend(loc="upper right", fontsize=14) # 凡例を表示
plt.tick_params(labelsize=14)
print("Figure 8: PSD of broadband data at 100-s period")
plt.show()
In [359]:
#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",fontsize=16)
# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(-120, -180)
plt.legend(loc="upper right", fontsize=14) # 凡例を表示
plt.tick_params(labelsize=14)
print("Figure 9: PSD of broadband data at 50-s period")
plt.show()
In [360]:
#plt.rcParams['figure.figsize'] = 11, 6
plt.rcParams['figure.figsize'] = 15, 5
plt.xlabel('Station',fontsize=16)
#ew_20db = pd.DataFrame({'x': [-9999, 99999],
# 'val': [20, 20]})
#plt.plot(ew_20db.x, ew_20db.val, label = 'Erhard Wielandt 20dB', color='black', linestyle='dashed')
sta_list = bb_100s.sta.as_matrix()
nz_val_med = bb_100s.nz_val_med.as_matrix()
ez_val_med = bb_100s.ez_val_med.as_matrix()
en_val_med = bb_100s.en_val_med.as_matrix()
x_position = np.arange(len(sta_list))
#print(x_position)
plt.bar(x_position - 0.3, nz_val_med, width=0.3, label='diff N-Z')
plt.bar(x_position, ez_val_med, width=0.3, label='diff E-Z')
plt.bar(x_position + 0.3, np.abs(en_val_med), width=0.3, label='abs. diff E-N')
#plt.bar(x_position, y_HNN, width=0.4, label='HNN')
#plt.bar(x_position + 0.4, y_HNE, width=0.4, label='HNE')
plt.legend(fontsize=16)
#plt.xticks(x_position + 0.2, x)
plt.xticks(x_position + 0.1, sta_list)
plt.ylabel('Relative PSD at 100-s period', fontsize=16)
#plt.ylabel('PSD at 200-s period', fontsize=16)
plt.title("Byerly Vault Relative PSDs at 100 s period",fontsize=16)
# HN
#plt.ylim(-85, -120)
# HH
plt.ylim(0, 30)
plt.legend(loc="upper right", fontsize=14) # 凡例を表示
plt.tick_params(labelsize=14)
print("Figure 10: Relative PSD of broadband data at 100-s period")
plt.show()
