美文网首页
获取波形任意时间上的振幅

获取波形任意时间上的振幅

作者: SeisBird | 来源:发表于2017-12-05 10:08 被阅读0次

    SAC提供的命令可以帮助用户实现地震数据的预处理,但无法实现所有的功能。比如我利用ppk进行震相走时标定之后,由于工作的需要,可能还需要知道该点的振幅大小。这就需要自己能够在别的程序中读写SAC文件,下面是利用Python的实现(你可以随意复制和修改,但我不承担可能造成的风险的责任,所以你最好需要读懂程序,读懂的前提是需要对SAC头段有一定的了解):

    #!/usr/bin/env python
    # This code is mainly used to output the amplitude of arbitrary time in SAC waveform files
    # Usage:
    #    python readsac.py dirname
    # Example:
    #    python readsac.py ~/project/sacfiles/ 
    # The source code is made by SeisBird, 2017/12/03. You can copy and modify in your ways. 
    # What you should know is that I am not responsibility for any influence on you for your modification.
    
    import os
    import re
    import sys
    import glob
    import numpy as np
    import matplotlib.pyplot as plt
    
    if len(sys.argv)!=2:
        sys.exit("Usage: python %s dirname\n" %sys.argv[0])
    
    dir=sys.argv[1]
    os.chdir(dir)
    for fname in glob.glob("*.sac"):
    
        f=open(fname,'r')
        lines=f.readlines()
        o=float(re.split('[\s]+',lines[1])[3])
        delta=float(re.split('[\s]+',lines[0])[1])
        npts=float(re.split('[\s]+',lines[15])[5])
    
        o_range=np.arange(o,o+delta*npts,delta)
        o_range_list=o_range.tolist()
    
        line=30
        data=[]
    
        while(line<len(lines)):
            data=data+re.split('[\s]+',lines[line])[1:6]
            line+=1
    
    
        if len(data)==npts:
            pass
        else:
            data.remove('')
    
        time_ampl=dict(zip(o_range_list,data))
    #   print(time_ampl)
        print(fname)
        while True:
            time=float(input("Please input time: "))
            j=0
        
            if time<0:
                break
            elif time in o_range_list:
                print("amplitude in %f is %s" %(time,time_ampl[time]))
            else:
                while(j<len(o_range_list)):
         
                    if(time>o_range_list[j] and time<o_range_list[j+1]):
                        amplitude=(float(time_ampl[o_range_list[j+1]])-float(time_ampl[o_range_list[j]]))/(o_range_list[j+1]-o_range_list[j])*(time-o_range_list[j])+float(time_ampl[o_range_list[j]])
                        print("amplitude in %f is %f" %(time,amplitude))
                        j+=1
                    else:
                        j+=1
    
        plt.plot(o_range_list,data)
        plt.show()
    
        print("The Next file:")
    

    运行之后,程序读取文件,输入自变量时间,得到振幅数据,如果你想结束进入到下一个文件,只需要输入的时间为负数即可,此时会产生波形的图像,可以检验。
    P.S. 这个程序写的匆忙,简单粗暴,着急用,而且我猜应该没多少人感兴趣,写的随便,就当记录了,所以可能会有一些bug和需要优化的地方,欢迎提出,联系我:seisbird@gmail.com

    修改历史

    1. 原稿 2017年12月4日

    相关文章

      网友评论

          本文标题:获取波形任意时间上的振幅

          本文链接:https://www.haomeiwen.com/subject/lmjqixtx.html