自己開發的一個SEED格式地震數據轉換為SAC格式數據，並完成世界時整天波形合並的Python腳本

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