Tidal signal measurement¶
This will download seismic waveforms & measure sensor orientation for a target station. ObsPy, numpy, scipy, matplotlib, pandas are required.
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__)
Import SciPy, NumPy, matplotlib, pandas module¶
In [2]:
import numpy as np
import scipy as sp
import matplotlib as mpl
import pandas as pd
print("# numpy version = ",np.__version__)
print("# scipy version = ",sp.__version__)
print("# matplotlib version = ",mpl.__version__)
print("# pandas version = ",pd.__version__)
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.dates as mdates
from matplotlib.lines import Line2D
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
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
function: get_seedid¶
to extract seed_id from obspy stream.
In [4]:
def get_seedid (tr):
seedid=tr.stats.network+"."+tr.stats.station+"."+tr.stats.location+"."+tr.stats.channel
return seedid
function: st_remove_resp¶
to correct instrument response.
In [5]:
def st_remove_resp (st, deciopt, decifactor, pre_filt, output):
st.detrend() # remove linear trend
st.detrend("demean") # demean
st.taper(0.05) # cosin taper
if deciopt == 1:
# decimate to 100Hz
if decifactor == 100:
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
elif decifactor == 300:
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
st.decimate(3, strict_length=False)
elif decifactor == 1000:
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
st.decimate(5, strict_length=False)
st.decimate(2, strict_length=False)
else:
st.decimate(factor=decifactor, strict_length=False)
st = st.remove_response(pre_filt=pre_filt,output=output,water_level=None) # get velocity data (m/s)
return st
function: dB¶
compute dB
In [6]:
def dB(x, out=None):
if out is None:
return 10 * np.log10(x)
else:
np.log10(x, out)
np.multiply(out, 10, out)
def dBpower(x, out=None):
if out is None:
return 10 * np.log10(x)
else:
np.log10(x, out)
np.multiply(out, 20, out)
function: nextpow2¶
In [7]:
def nextpow2(i):
n = 1
while n < i: n *= 2
return n
Theoretical tidal periods¶
In [8]:
tide_theo_hr = {"M2": 12.42060, "S2": 12.00000, "N2": 12.65835, "K2": 11.96724, "K1": 23.93447, "O1": 25.81934, "P1": 24.06589, "Q1":26.86836 }
print(tide_theo_hr)
Set time window for downloading seismic data¶
This example extracts ~16-month-long data
In [9]:
starttime = UTCDateTime("2007-06-30T12:00:00")
endtime = UTCDateTime("2008-11-08T00:00:00")
#print("# starttime = ",starttime)
#print("# endtime = ", endtime)
elapsed_day = (endtime - starttime) / (60*60*24)
#print("# elapsed_day = ",elapsed_day)
Set SNCL parameters¶
Which SNCL (Station, Network, Component, Location) we will use? sta1,net1,com1,loc1 for the target data.
In [10]:
# targert
sta1 = "YBH" # station
net1 = "BK" # network
com1 = "LHZ" # component "HH?" to get all 3-com data from broadband sensor
loc1 = "--" # location
deciopt_1 = 1
decifactor_1 = 300 # 1sps -> 0.03sps
client1 = Client("NCEDC") # data from NCEDC
Directory for waveform plot¶
This example will create directory "YBH_tidal_analysis" where all plots will be saved.
In [11]:
# name for output directory
pwd_dir = os.getcwd()
#sacdir= pwd_dir +"/"+ evmseedid +"_M"+(str)(evmag)+"_"+(str)(eventid_ncedc)+"_fl"+(str)(f2)+"_fh"+str(f3)+"_dist"+(str)(distkm_from_eq)+"km"
#plot_dir= pwd_dir +"/"+ evmseedid +"_M"+(str)(evmag)+"_"+(str)(evid)
#+"_fl"+(str)(f2)+"_fh"+str(f3)+"_dist"+(str)(distdeg_from_eq)+"deg"
plot_dir= pwd_dir +"/"+sta1+"_tidal_analysis"
print("# plot_dir = ",plot_dir)
# create output directory
if not os.path.exists(plot_dir):
os.makedirs(plot_dir)
Set frequency range for sensor orientation¶
set the frequency range (fl and fh) for correcting the instrument response. This example uses 100-day to 1.0 Hz band. fl_wave & fh_wave are a filter band (0.01-1mHz) for plotting gravitational tidal force.
In [12]:
#fl2 < fl < fh < fh2
fl2 = 1.0 / (24*60*60*120) # 1/120 day
fl = 1.0 / (24*60*60*100) # 1/100 day
fh = 1.0
fh2 = 2.0
# waveform plot
fl_wave = 0.00001
fh_wave = 0.001
fl_wave_mHz = fl_wave * 1000 # convert Hz -> mHz
fh_wave_mHz = fh_wave * 1000 # convert Hz -> mHz
#pre_filt = [0.015, 0.02, 45, 50]
pre_filt = [fl2, fl, fh, fh2]
Set parameters for Welch spectral estimation¶
This example uses 400-day-long data with 80% overlapping for Welch spetracal estimation
In [13]:
w_length_day = 400 # days
overlap = 0.80 # 80%
overlap_percent = overlap*100
In [14]:
welch_para = "Total "+str(elapsed_day )+"-day-long data. Amplitude spectra: "+str(w_length_day)+"-day time window with "+str(overlap_percent)+"% overlapping"
print("# welch_para = ", welch_para)
Downloading seismic data¶
use get_waveforms to download data and do st.plot() for plotting. Also use get_stsations to obtation station inventory file.
In [15]:
st1 = client1.get_waveforms(network=net1, station=sta1, location=loc1, channel=com1,
starttime=starttime, endtime=endtime,
attach_response=True)
inv1 = client1.get_stations(network=net1, station=sta1, location=loc1, channel=com1,
starttime=starttime, endtime=endtime,
level="response")
_plot = st1.plot(size=(1000,250))
#_plot = st1.plot()
Removing instrument response¶
use st_remove_resp function uses obspy remove_response to remove instrument response. Example will provide ground acceleration data (m/s**2)
In [16]:
st1_gm = st_remove_resp(st1.copy(), deciopt_1, decifactor_1, pre_filt, "ACC")
_plot = st1_gm[0].plot(size=(1000,250))
Theoretical gravity data¶
This data is calculated by SPOTL software (https://igppweb.ucsd.edu/~agnew/Spotl/spotlmain.html). Channel name for gravity data is RHG.
In [17]:
grav_fi = "http://ncedc.org/ftp/outgoing/taira/WQC/YBH.grav.out"
grav =pd.read_csv(grav_fi,
names=["microgal"],
header=None, skiprows=0)
#type(grav['microgal'])
grav_out = (grav['microgal']*1.0).to_numpy()
st1_grav = st1_gm.copy()
st1_grav[0].data = grav_out
st1_grav[0].stats.channel = "RHG"
st1_grav[0].stats.delta = 300
st1_grav[0].stats.sampling_rate = 1/st1_grav[0].stats.delta
st1_grav[0].stats.starttime = starttime
In [18]:
st1_merge = st1_gm.copy() + st1_grav.copy()
#print(st1_merge)
Checking if each segment has enough data points¶
In [19]:
for tr in st1_merge:
npts_out=tr.stats.npts
sta_out = tr.stats.station
sampling_rate_out = tr.stats.sampling_rate
print("# sta_out = ", sta_out, " npts_out = ", npts_out)
w_length = int(w_length_day * 60 * 60 * 24 * sampling_rate_out)
print("# sampling_rate_out = ",sampling_rate_out," w_lengh = ", w_length)
if npts_out < w_length:
print("# trace removed low length. npts_out = ", npts_out)
st1_merge.remove(tr)
In [20]:
st1_merge.detrend() # remove linear trend
st1_merge.detrend("demean") # demean
st1_merge.taper(0.05) # cosin taper
Out[20]:
Checking waveforms¶
In [21]:
st1_select = st1_merge.copy()
t1 = UTCDateTime("2008-04-01T00:00:00")
t2 = UTCDateTime("2008-05-01T00:00:00")
st1_select.detrend() # remove linear trend
st1_select.detrend("demean") # demean
st1_select.taper(0.05) # cosin taper
#print("# fl_wave = ",fl_wave)
st1_select.filter("bandpass", freqmin = fl_wave , freqmax = fh_wave, corners=6, zerophase=True)
for tr in st1_select:
#tr.plot()
tr.plot(starttime=t1, endtime=t2, size=(1000,250))
Computing Welch spectra¶
In [22]:
i = 0
for tr in st1_merge:
print("# tr = ",tr)
delta = tr.stats.delta
#
#w_length = 1728000
#w_length = 172800*2
#w_length = 172800*3
sta_out = tr.stats.station
sampling_rate_out = tr.stats.sampling_rate
npts_out = tr.stats.npts
#print("# sta_out = ", sta_out, " npts_out = ", npts_out)
w_length = int(w_length_day * 60 * 60 * 24 * sampling_rate_out)
#print("# sampling_rate_out = ",sampling_rate_out," w_lengh = ", w_length)
nperseg = w_length
noverlap = int(w_length * overlap)
#nperseg = tr.stats.npts
#print("# nperseg = ", nperseg)
nfft = nextpow2(nperseg)
#print("# nfft = ", nfft)
#print("# delta = ",delta)
#for i, tr in enumerate(st):
#print(" type = ", type(tr))
#print("# i = ",i)
#fre1_sp, p11_sp = signal.welch(tr, 1/delta, nperseg=nperseg, noverlap=noverlap)
#freq, p11 = signal.welch(tr, 1/delta, nperseg=nperseg, noverlap=noverlap)
#freq, p11 = signal.welch(tr, 1/delta, nperseg=nperseg, noverlap=noverlap, nfft=57600)
#freq, p11 = signal.welch(tr, 1/delta, nperseg=nperseg, noverlap=noverlap, nfft=147744)
freq, p11 = signal.welch(tr, 1/delta, nperseg=nperseg, noverlap=noverlap, nfft=nfft)
tr.fdata = p11
tr.freq = freq
tr.p11 = p11
tr.p11_freq = freq
#print("# tr.fdata = ",tr.fdata)
#print("# i = , i")
#print(st1[0].fdata)
Plotting waveforms and amplitude spectra¶
In [23]:
offset = 160 # convert m/s^2 to microGal. 100 m/s^ = 1 Gal. 20xlog10(10^8)
seis_seedid=st1_merge[0].stats.network+"."+st1_merge[0].stats.station+"."+st1_merge[0].stats.location+"."+st1_merge[0].stats.channel
#grav_seedid=st1_merge[1].stats.network+"."+st1_merge[1].stats.station+"."+st1_merge[1].stats.location+"."+st1_merge[1].stats.channel[:-1]
#plot_xc.BRIB.73292045.BRIB_ref.pdf
#plot_fi = plot_dir+"/plot_xc_"+target_seedid+"_ref"+ref_seedid+"_"+event_para+".pdf"
#plot_fi = "test.pdf"
plot_fi = plot_dir+"/plot_tidal_signal_"+seis_seedid+".pdf"
gs = gridspec.GridSpec(3,2)
#gs = gridspec.GridSpec(6,2)
#fig = plt.figure()
fig = plt.figure(figsize=(16, 16))
plt.subplots_adjust(wspace=0.45, hspace=0.4)
data_window_day = 35 # day
xdata_point_min = 100000 # starting data point
data_window_sec = data_window_day * (60*60*24)
data_point = data_window_sec * sampling_rate_out
#data_all = np.concatenate((st[0].data, st[1].data) )
#xdata_point_max = 110000
xdata_point_max = xdata_point_min + int(data_point)
sampling_rate_out = st1_select[0].stats.sampling_rate
#print("# sampling_rate_out = ", sampling_rate_out)
data_window_sec = (xdata_point_max - xdata_point_min) / sampling_rate_out
#print("# data_window_sec = ", data_window_sec)
data_window_day = data_window_sec / (60*60*24)
#print("# data_window_day = ", data_window_day)
data_window_day_out = "{:.2f}".format(data_window_day)
t=fig.text(0.13, 0.86, "Selected "+str(data_window_day_out)+"-day data")
t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))
t=fig.text(0.13, 0.84, str(fl_wave_mHz)+"-"+str(fh_wave_mHz)+" mHz BP filter")
t.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='white'))
#plt.subplot(gs[0, 0])
ax1 = plt.subplot(gs[0, :])
#plt.title(event_para+" \n Origin Time:"+str(origin_time)+" USGS event-id:"+evid \
# +" \n "+welch_para)
plt.title("Seismic data: "+seis_seedid+"\n Total "+str(elapsed_day )+"-day-long data starting from "+str(starttime)+"\n Amplitude spectra: "+str(w_length_day)+"-day time window with "+str(overlap_percent)+"% overlapping")
#plt.title(welch_para)
sncl1 = get_seedid(st1_select[0])
sncl2 = get_seedid(st1_select[1])
data_all = st1_select[0].data
wf_amp_max = max(data_all[xdata_point_min:xdata_point_max].min(), data_all[xdata_point_min:xdata_point_max].max(), key=abs)
if wf_amp_max < 0:
wf_amp_max = wf_amp_max * -1
#print("# wf_amp_max = ", wf_amp_max)
wf_amp_max2 = wf_amp_max * 1.35
wf_amp_min2 = wf_amp_max2 * -1
wf_amp_max2_microgal = wf_amp_max2 * 100000000 # to gonvert m/s^2 -> microGal
wf_amp_min2_microgal = wf_amp_max2_microgal * -1
#print("# wf_amp_max2_microgal = ", wf_amp_max2_microgal)
ax1.plot_date(st1_select[0].times("matplotlib"), st1_select[0].data*1, fmt='-', label=sncl1, color="red", linewidth=1.0, linestyle='solid')
#ax1.set_ylabel("Amplitude (m/s2)")
ax1.set_ylabel("Amplitude (m/s$^2$)")
#ax1.set_xlim(st1_select[0].times("matplotlib")[0], st1_select[0].times("matplotlib")[-1])
ax1.set_xlim(st1_select[0].times("matplotlib")[100000], st1_select[0].times("matplotlib")[110000])
ax1.set_ylim(wf_amp_min2 , wf_amp_max2 )
#myFmt = mdates.DateFormatter("%D/%H:%M")
#myFmt = mdates.DateFormatter("%Y-%m-%d\n%H:%M:%S")
myFmt = mdates.DateFormatter("%Y\n%m/%d")
#print("# myFmt = ", myFmt)
#plt.gca().xaxis.set_major_formatter(myFmt)
#ax1.xaxis.set_major_formatter(myFmt)
plt.grid()
ax2 = ax1.twinx()
ax2.plot_date(st1_select[1].times("matplotlib"), st1_select[1].data*1, fmt='-', label=sncl2, color="blue", linewidth=1.0, linestyle='dashed')
#ax2.set_ylim([200, -150])
ax2.set_ylim(wf_amp_max2_microgal , wf_amp_min2_microgal )
#plt.title(event_para+" \n Origin Time:"+str(origin_time)+" USGS event-id:"+evid+" \n M"+str(evmag)+" "+event_region+" \n Target station: "+target_seedid+" Ref. station: "+ref_seedid+" \n Diff. dist.: "+diff_eq_dist_out+" Diff. az.: "+diff_eq_baz_out)
ax2.xaxis.set_major_formatter(myFmt)
ax2.set_ylabel("Gravity ($\mu$Gal)")
#plt.xlabel("Time")
#plt.legend(loc = "upper right")
#plt.show()
colors = ['red', 'blue']
styles = ['solid', 'dashed']
#lines = [Line2D([0], [0], color=c, linewidth=1, linestyle="solid") for c in colors]
#for i, color in enumerate(colors):
#lines = [Line2D([0], [0], color="black", linewidth=1.5, linestyle=s) for s in styles]
lines = [Line2D([0], [0], color=colors[i], linewidth=1.5, linestyle=s) for i, s in enumerate(styles) ]
#labels = ['Residual energy', 'Cross-correlation value' ]
#labels = [sncl1, sncl2]
labels = ['Observed seismic data', 'Theoretical gravity tide']
plt.legend(lines, labels, loc = "upper right")
### PSD ###
ax1 = plt.subplot(gs[1,:])
for k in tide_theo_hr.keys():
#print("# k = ", k)
val=tide_theo_hr[k]
#print("# val = ", val)
plt.vlines(x=val , ymin=-999, ymax=999, color="black", linewidth = 0.75, linestyle = "--")
# seismic
ax1.plot( (1.0/st1_merge[0].freq[1:])/(60*60), dB(st1_merge[0].fdata[1:]), linewidth = 1.0, color="red" )
ax1.set_ylim(-150,-50)
#ax1.set_ylabel("Amplitude (m$^2$/s$^4$/Hz)")
ax1.set_ylabel("Power [10log$_{10}$(m$^2$/s$^4$/Hz)](dB)")
ax1.set_xlabel("Period (hour)")
plt.grid()
ax2 = ax1.twinx()
# gravity
ax2.plot( (1.0/st1_merge[1].freq[1:])/(60*60), (dB(st1_merge[1].fdata[1:])), linewidth = 1.0, linestyle="dashed", color="blue" )
ax2.set_ylim(-150+offset,-50+offset)
#ax2.set_ylabel("Gravity ($\mu$Gal$^2$)")
ax2.set_ylabel("Power [10log$_{10}$($\mu$Gal$^2$/Hz)](dB)")
ax2.set_xlim([10, 35])
#ax2.set_xlim([10, 13])
#ax2.set_xlim([22, 28])
#plt.xlabel("Period (hour)")
colors = ['red', 'blue']
styles = ['solid', 'dashed']
#lines = [Line2D([0], [0], color=c, linewidth=1, linestyle="solid") for c in colors]
#for i, color in enumerate(colors):
#lines = [Line2D([0], [0], color="black", linewidth=1.5, linestyle=s) for s in styles]
lines = [Line2D([0], [0], color=colors[i], linewidth=1.5, linestyle=s) for i, s in enumerate(styles) ]
#labels = ['Residual energy', 'Cross-correlation value' ]
#labels = [sncl1, sncl2]
labels = ['Observed seismic data', 'Theoretical gravity tide']
plt.legend(lines, labels, loc = "upper right")
### PSD ###
### PSD around 12-hour ###
ax1 = plt.subplot(gs[2,0])
for k in tide_theo_hr.keys():
#print("# k = ", k)
val=tide_theo_hr[k]
#print("# val = ", val)
plt.vlines(x=val , ymin=-999, ymax=999, color="black", linewidth = 0.75, linestyle = "--")
# seismic
ax1.plot( (1.0/st1_merge[0].freq[1:])/(60*60), dB(st1_merge[0].fdata[1:]), linewidth = 1.0, color="red" )
ax1.set_ylim(-150,-50)
ax1.set_ylabel("Power [10log$_{10}$(m$^2$/s$^4$/Hz)](dB)")
ax1.set_xlabel("Period (hour)")
plt.grid()
ax2 = ax1.twinx()
# gravity
ax2.plot( (1.0/st1_merge[1].freq[1:])/(60*60), (dB(st1_merge[1].fdata[1:])), linewidth = 1.0, linestyle="dashed", color="blue" )
ax2.set_ylim(-150+offset,-50+offset)
ax2.set_ylabel("Power [10log$_{10}$($\mu$Gal$^2$/Hz)](dB)")
ax2.set_xlim([11, 14.5])
colors = ['red', 'blue']
styles = ['solid', 'dashed']
lines = [Line2D([0], [0], color=colors[i], linewidth=1.5, linestyle=s) for i, s in enumerate(styles) ]
labels = ['Seismic data', 'Gravity data']
plt.legend(lines, labels, loc = "upper right")
### PSD around 12-hour ###
### PSD around 24-hour ###
ax1 = plt.subplot(gs[2,1])
for k in tide_theo_hr.keys():
#print("# k = ", k)
val=tide_theo_hr[k]
#print("# val = ", val)
plt.vlines(x=val , ymin=-999, ymax=999, color="black", linewidth = 0.75, linestyle = "--")
# seismic
ax1.plot( (1.0/st1_merge[0].freq[1:])/(60*60), dB(st1_merge[0].fdata[1:]), linewidth = 1.0, color="red" )
ax1.set_ylim(-150,-50)
ax1.set_ylabel("Power [10log$_{10}$(m$^2$/s$^4$/Hz)](dB)")
ax1.set_xlabel("Period (hour)")
plt.grid()
ax2 = ax1.twinx()
# gravity
ax2.plot( (1.0/st1_merge[1].freq[1:])/(60*60), (dB(st1_merge[1].fdata[1:])), linewidth = 1.0, linestyle="dashed", color="blue" )
ax2.set_ylim(-150+offset,-50+offset)
ax2.set_ylabel("Power [10log$_{10}$($\mu$Gal$^2$/Hz)](dB)")
ax2.set_xlim([20, 34.5])
colors = ['red', 'blue']
styles = ['solid', 'dashed']
lines = [Line2D([0], [0], color=colors[i], linewidth=1.5, linestyle=s) for i, s in enumerate(styles) ]
labels = ['Seismic data', 'Gravity data']
plt.legend(lines, labels, loc = "upper right")
### PSD around 24-hour ###
#plt.savefig("test.pdf")
plt.savefig(plot_fi)
In [ ]:
In [ ]:
In [ ]:
In [ ]: