import obspy
import warnings
import os
import time
from numba import jit
warnings.filterwarnings("ignore")
delta = 100 # 仪器采样率
# 国产SEED格式文件为BJT，文件开始时间为文件名前一天的16时，因此合并后的第1个完整的世界时以天为单位的数据是目录中第2个
# 文件名中日期，要注意要想转换合并一年的完整数据，要在目录中加入上一年最后一个BJT的SEED文件，和下一年中第1个BJT文件。
data_path = r'D:\QXZ\2020' # 需要转换格式，合并数据的文件目录，此目录中存放的是SEED格式原始数据文件
output_dir = r'D:\JL\QXZ\2020' # 合并后文件的保存路径
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# 判断两个数据文件的日期是否连续
@jit(forceobj=True) # Set "nopython" mode for best performance
def merge_st(st):
print('st.merge starting....')
st.merge(fill_value=0) # 国产SEED格式文件有问题，直接转存为SAC时会分隔成很多小的文件，主要是数据重叠或断记产生的，此命令完成
# 各通道数据整合，并将断记点自动填充为0，解决了转存后会产生很多小文件的问题
print('st.merge completed!')
def is_continuoues(first_event, second_event):
try:
st1 = obspy.read(data_path + '\\' + first_event)
first_date = str(st1[0].stats.starttime).split('T')[0]
st2 = obspy.read(data_path + '\\' + second_event)
second_date = str(st2[0].stats.starttime).split('T')[0]
except:
return(False, 1, 2)
first_julday = obspy.UTCDateTime(first_date + 'T00:00:00.000000Z').timestamp
second_julday = obspy.UTCDateTime(second_date + 'T00:00:00.000000Z').timestamp
if second_julday - first_julday == 86400.0:
return (True, st1, st2)
else:
return (False, st1, st2)
rest_event = os.listdir(data_path)
# 遍历数据目录
for event in os.listdir(data_path):
if len(rest_event) == 1:
print('完成全部数据合并，程序结束')
break
first_event = event
rest_event.pop(0)
second_event = rest_event[0]
is_cont, st1, st2 = is_continuoues(first_event, second_event)
if is_cont:
print(f'开始合并{second_event}............')
t_start = time.time()
start_time = str(st2[0].stats.starttime).split('T')[0]+'T00:00:00.000000Z'
end_time = str(st2[0].stats.starttime).split('T')[0]+'T23:59:59.990000Z'
start_point = 28800 * delta # 截取一天波形起始点
st = st1 # 读取第1个数据文件
st += st2 # 将第2个数据文件合并到第1个数据文件中，注意现在的st变量仍是SEED格式。
st = merge_st(st)
st[0].data = st[0].data[start_point:(start_point+8640000)]
st[0].stats.starttime = obspy.UTCDateTime(start_time)
st[0].stats.npts = 8640000
st[1].data = st[1].data[start_point:(start_point+8640000)]
st[1].stats.starttime = obspy.UTCDateTime(start_time)
st[1].stats.npts = 8640000
st[2].data = st[2].data[start_point:(start_point+8640000)]
st[2].stats.starttime = obspy.UTCDateTime(start_time)
st[2].stats.npts = 8640000
net = st[0].stats.network
sta = st[0].stats.station
year = start_time[0:4]
jday = obspy.UTCDateTime(start_time).julday
if jday < 10:
jday = '00' + str(jday)
elif jday < 100:
jday = '0' + str(jday)
jday = str(jday)
loc = st[0].stats.location
chan = []
chan.append(st[0].stats.channel)
chan.append(st[1].stats.channel)
chan.append(st[2].stats.channel)
file_name = []
file_name.append(sta+'.'+net+'.'+loc+'.'+str(chan[0])+'.'+year+'.'+jday)
file_name.append(sta+'.'+net+'.'+loc+'.'+str(chan[1])+'.'+year+'.'+jday)
file_name.append(sta+'.'+net+'.'+loc+'.'+str(chan[2])+'.'+year+'.'+jday)
st[0].write(output_dir+'\\'+file_name[0], format='SAC')
st[1].write(output_dir+'\\'+file_name[1], format='SAC')
st[2].write(output_dir+'\\'+file_name[2], format='SAC')
t_end = time.time()
print(f'已经生成文件{file_name[0]}!用时{t_end-t_start}秒')
else:
print(f'数据不连续！{second_event}合并失败！')
continue
