import os
import numpy as np
import pandas as pd
import scipy as sp
import matplotlib.pyplot as plt
import pickle
from glob import glob
from obspy import read
from scipy.signal import (
butter,
detrend,
filtfilt,
tukey,
)
!ls
2021-03-04T19:28:32.820000_M_8.10 2022-09-19T18:05:06.320000_M_7.60 2021-07-29T06:15:48.030000_M_8.20 plot_DAS_teleseismic_earthquake.ipynb 2021-08-14T12:29:08.420000_M_7.20 teleseismic_earthquakes.csv 2021-09-08T01:47:47.930000_M_7.00
# load station information
sta_df = pd.read_csv('../das_and_station_location/stations_location.csv')
print(sta_df)
net sta lat lon ele 0 CI SRT 35.69235 -117.75051 667.0 1 CI LRL 35.47954 -117.68212 1340.0 2 CI MLAC 37.63019 -118.83605 2162.0 3 NP 5419 35.64906 -117.66187 689.0
# load das channel information
das_df = pd.read_csv('../das_and_station_location/Ridgecrest_North_DAS_ChannelLocation.txt')
das_good_df=das_df[das_df['status']=='good'].reset_index(drop=True)
print(das_good_df)
channel latitude longitude elevation x status 0 146 35.651884 -117.669510 NaN 0.000000 good 1 147 35.651884 -117.669599 NaN 0.008037 good 2 148 35.651885 -117.669688 NaN 0.016075 good 3 149 35.651885 -117.669777 NaN 0.024112 good 4 150 35.651885 -117.669866 NaN 0.032150 good ... ... ... ... ... ... ... 1021 1244 35.652205 -117.760332 NaN 8.206211 good 1022 1245 35.652205 -117.760421 NaN 8.214248 good 1023 1246 35.652205 -117.760510 NaN 8.222286 good 1024 1247 35.652206 -117.760599 NaN 8.230323 good 1025 1248 35.652206 -117.760688 NaN 8.238360 good [1026 rows x 6 columns]
# load the event information
catalog = pd.read_csv('teleseismic_earthquakes.csv', sep=' ')
event=catalog.iloc[0]
print(event)
Magnitude 8.2 Time 2021-07-29T06:15:48.030 Latitude 55.4516 Longitude -157.964 Depth 28780.0 Name: 0, dtype: object
event_path=glob(f"{event['Time']}*")[0]
print(event_path)
print(glob(f"{event_path}/*"))
2021-07-29T06:15:48.030000_M_8.20 ['2021-07-29T06:15:48.030000_M_8.20/CI_SRT_VEL.mseed', '2021-07-29T06:15:48.030000_M_8.20/Ridgecrest_North_das_info.pkl', '2021-07-29T06:15:48.030000_M_8.20/Ridgecrest_North_das_strain.npy']
# load sta waveforms and plot
st = read(glob(f"{event_path}/*.mseed")[0])
st.plot()
# load das waveforms information
info_data_das = pickle.load(open(glob(f"{event_path}/*_das_info.pkl")[0], 'rb'))
for key in info_data_das.keys():
print(key, ':', info_data_das[key])
dx : -1.0 dt : 0.02 fs : 50.0 nt : 270000 nx : 1026 begTime : 2021-07-29 06:15:48.034400 endTime : 2021-07-29 07:45:48.014400 raw2strain_or_strainrate_done : True unit : strain: 1(=m/m)
fs = 50 Hz is the sampling rate
nx = 1026 is the total channel number, same as it in 'das_good_df', channel location can be found in 'das_good_df'.
# load das waveforms data
data_das = np.load(glob(f"{event_path}/*_das_strain.npy")[0])
# preprocess das data
data = detrend(data_das, axis=-1, type="constant")
data = detrend(data, axis=-1, type="linear")
window = tukey(data.shape[-1], alpha=0.05)
data = data * window
# bp filter
fs = info_data_das['fs']
freqmin = 0.01 #Hz
freqmax = 0.1 #Hz
passband = [freqmin * 2 / fs, freqmax * 2 / fs]
b, a = butter(2, passband, "bandpass")
data_bp = filtfilt(b, a, data)
# bp the STA data
st_bp = (
st[0]
.copy()
.filter('bandpass', freqmin=freqmin, freqmax=freqmax, zerophase=True, corners=2)
)
# plot das data
plt.figure(figsize=(8, 4))
vmax = np.percentile(np.abs(data_bp), 99)
plt.imshow(
data_bp,
extent=(0, 5400, data_bp.shape[0], 0),
aspect='auto',
cmap='seismic',
vmax=vmax,
vmin=-vmax,
)
vmax2 = np.max(np.abs(st_bp.data))
plt.plot(
np.linspace(0, 5400, st_bp.stats.npts),
data_bp.shape[0] / 2 + st_bp.data / vmax2 * data_bp.shape[0] / 3,
'k',
label='Station',
)
plt.legend()
plt.xlim(0, 1500)
plt.xlabel(f'Time after origin={event["Time"]} (s)')
plt.ylabel('Channel')
plt.show()