Python-SAC-IO

Posted on | In

简介

在Python中使用包实现SAC文件读取。有如下三个SAC读取相关的库：
pysmo
sacpy
pyrocko

注意：有些库不支持较高的Python版本。例如pyrocko只支持到3.6.*版本。使用conda创建环境较为方便，建议使用3.6版本的
python，在改版本下均可安装以下库及Obspy。

运行如果出现报错：

'version_info' object has no attribute '__version__'

试一下运行：

pip install pyparsing==2.4.7

pysmo

安装

pip install pysmo
pip install pysmo --pre
pip install git+https://github.com/pysmo/pysmo

使用

官方文档：https://pysmo.readthedocs.io/en/latest/sac.html

读取文件

1
2
3
4
from pysmo import SacIO
seismogram = SacIO.from_file('filename.sac')
datas = seismogram.data
print(datas)

读取采样率

1
2
delta = seismogram.delta
print(delta)

修改采样率

1
seismogram.delta = 0.05

读取IRIS的地震数据

1
2
3
4
5
6
7
8
9
10
11
12
13
>>> from pysmo import SacIO
>>> my_sac = SacIO.from_iris(
>>>             net="C1",
>>>             sta="VA01",
>>>             cha="BHZ",
>>>             loc="--",
>>>             start="2021-03-22T13:00:00",
>>>             duration=1 * 60 * 60,
>>>             scale="AUTO",
>>>             demean="true",
>>>             force_single_result=True)
>>> my_sac.npts
144001

sacpy

安装

使用pip即可

pip install sacpy

pyrocko

Obspy

Obspy也可以直接读取sac文件。如下：

1
2
3
4
5
6
7
8
9
10
import obspy
import numpy as np
import matplotlib.pyplot as plt

st = obspy.read("filename.sac")
tr=st[0]
print(st)
print(type(st))
print(st[0].data)
print(st[0].stats.delta)

注意： 使用obspy读取得到的为stream格式，该格式是一个列表。需要使用st[0]表示sac的图像数据。

获取基本信息

st[0].stats

1
2
3
4
5
6
7
8
9
10
11
12
      network: NA
      station: BILL
     location: 
      channel: LHE
    starttime: 2000-01-14T23:37:10.800000Z
      endtime: 2000-01-15T00:30:29.800000Z
sampling_rate: 1.0
        delta: 1.0
         npts: 3200
        calib: 1.0
      _format: SAC
          sac: AttribDict({'delta': 1.0, 'depmin': -4299.8511, 'depmax': 4099.0977, 'b': 0.0, 'e': 3199.0, 'o': 0.0, 'internal0': 2.0, 'stla': 68.065102, 'stlo': 166.45239, 'stdp': 0.0, 'evla': 25.610001, 'evlo': 101.06, 'evdp': 33000.0, 'dist': 6370.9048, 'az': 23.937832, 'baz': 257.17334, 'gcarc': 57.340527, 'depmen': -2.0147384e-06, 'nzyear': 2000, 'nzjday': 14, 'nzhour': 23, 'nzmin': 37, 'nzsec': 10, 'nzmsec': 800, 'nvhdr': 6, 'norid': 0, 'nevid': 0, 'npts': 3200, 'iftype': 1, 'iztype': 9, 'leven': 1, 'lpspol': 0, 'lovrok': 1, 'lcalda': 1, 'unused23': 0, 'kstnm': 'BILL', 'kcmpnm': 'LHE', 'knetwk': 'NA', 'kevnm': ''})

获取采样率

st[0].stats.delta

去趋势

使用：

st=obspy.read(filepath)
st[0].detrend('polynomial', order=3) # 进行去趋势，选择为多项式，阶数为3

去除线性或者多阶趋势。
选项1可选：

  • simple
    Subtracts a linear function defined by first/last sample of the trace (uses obspy.signal.detrend.simple()).
  • linear
    Fitting a linear function to the trace with least squares and subtracting it (uses scipy.signal.detrend()).
  • constant or demean
    Mean of data is subtracted (uses scipy.signal.detrend()).
  • polynomial
    Subtracts a polynomial of a given order. (uses obspy.signal.detrend.polynomial()).
  • spline
    Subtracts a spline of a given order with a given number of samples between spline nodes.(uses obspy.signal.detrend.spline()).

示例代码

使用方法

在最后两行：
第一行为图所在目录（文件夹），文件为SAC格式。如果不是SAC格式，将被忽略。
第二行为：

  • path1 路径
  • 2 图类型type
    • type=1 速度图
    • type=2 加速度图
      均转换为位移图。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
import numpy
import sacpy
import pysmo
# import pyrocko
import obspy
import numpy as np
import matplotlib.pyplot as plt
from pysmo import SacIO
import os
from obspy.signal.detrend import polynomial  # 去除多阶趋势

def cusdetrend(filepath, tr, order=3):  # 去趋势
# 文件路径；由读取sac文件得到的trace，即seismogram[0]；执行多项式去趋势的阶数，默认为3阶
    seismogram = tr.copy()
    been_detrend = seismogram.detrend('polynomial', order=order) # 进行去趋势，选择为多项式，阶数为3
    return been_detrend


def integral(input_file, dt):  # 执行积分
# 输入文件为列表文件；dt为时间采样间隔
    length = len(input_file)
    output = [0 for i in range(0, length)]
    for i in range(0, length-1):
        if i == 0:
            output[i] = 0
        else:
            time = 0
            for j in range(0, i):
                time = time + (input_file[j]+input_file[j+1])/2*dt
            output[i] = time

    return output

def cus2disp(input_file, dt, type): # 转换为位移图， 
# 输入文件（列表）；采样时间间隔；输入图形类型（1为速度图，2为加速度图）
    if type == 1:
        outfile = integral(input_file, dt)
    elif type == 2:
        outfile = integral(integral(input_file, dt), dt)
    else:
        print('wrong input type!')
    return outfile

def savejpg(filepath, tr, imgtype):  # 绘制图形
# 文件路径；trace；图形类型，目前没有用
    # 调整画布大小，进行绘制
    t = np.arange(0, (tr.stats.npts / tr.stats.sampling_rate), tr.stats.delta)
    fig = plt.figure(1, figsize=(16, 9), dpi=200)
    plt.subplot(2, 1, 2)
    plt.plot(t, tr.data, 'k')
    plt.xlabel('T(s)')
    plt.ylabel('Displacement(nm)')
    plt.title('had detrend'+' '+tr.stats.station+tr.stats.channel)
    # 判断文件是否存在, 存在，则删除文件
    savepath=filepath+'detrend'+tr.stats.station+'_'+tr.stats.channel+'.jpg'
    if(os.path.exists(savepath)):
        os.remove(savepath)
    # 保存
    print('###############################')
    print('saving jpg file.')
    plt.savefig(filepath+'detrend'+tr.stats.station+'_'+tr.stats.channel+'.jpg')
    print(f'The file has been saved as {filepath}{tr.stats.station}_{tr.stats.channel}.jpg')
    print('###############################')
    # plt.show()
    plt.close(fig)
    plt.close('all')

    # 判断文件是否存在, 存在，则删除文件
    savepath=filepath+'detrend'+tr.stats.station+'_'+tr.stats.channel+'.sac'
    if(os.path.exists(savepath)):
        os.remove(savepath)
    print('###############################')
    print('saving sac file...')
    tr.write(savepath, format='SAC')
    print(f'The file has been saved as {filepath}{tr.stats.station}_{tr.stats.channel}.jpg')
    print('###############################')

def acc2dis_detrend(path, type):  # 主程序：转换为displacement并去趋势
# 文件路径；type=1 为velo，type=2为acc
    for file_name in os.listdir(path):
        if(file_name[-3:]!='SAC'):  # 避免读取其他文件
            continue
        filepath = path+'/'+file_name
        seismogram = obspy.read(filepath)
        trace = seismogram[0]
        # 进行转换,返回值为列表
        trace_data = cus2disp(trace.data, dt=trace.stats.delta, type=type)
        # 列表导入到stream数据中,只有这样才能使用obspy的去趋势函数。
        trace.data = numpy.array(trace_data[0:len(trace_data)-1])
        print(trace.data)
        # 进行去趋势
        detrend_result = cusdetrend(filepath, trace, order=3)  # 进行去趋势。
        # 保存
        if type == 1:
            imgtype = 'velocity(nm/s)'
        elif type == 2:
            imgtype = 'accelate(nm/s2)'
        savejpg(filepath, detrend_result, imgtype)

# 文件所在位置
path1 = 'G:/Pythonworks/Obspy_data_solve/seismogram/0_200acc'
# 执行
acc2dis_detrend(path1, 2)

归一化

安装 sklearn 库

conda 安装：

1
conda install -c anaconda scikit-learn