Python中调用SAC

简单示例

下面的脚本展示了如何在 Python 中调用 SAC。

下载地址： 0.simple-script.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import subprocess

os.putenv("SAC_DISPLAY_COPYRIGHT", '0')

s = "fg seismo \n"
s += "lh evla kstnm \n"
s += "q \n"
subprocess.Popen(['sac'], stdin=subprocess.PIPE).communicate(s.encode())

Python 中使用 subprocess 模块的 Popen 方法调用 SAC，通过 p.communicate() 将命令 s.encode() 传递给 SAC。

数据转换

首先要将 SEED 格式的数据转换成 SAC 格式。

• rdseed 一次只能处理一个 SEED 数据

• rdseed 的 -pdf 选项会提取出 SAC 波形数据和 PZ 格式的仪器响应文件

下载地址： 1.rdseed.py

警告

IRIS 自 2020 年 1 月开始不再提供 SEED 格式的数据下载支持，详情可阅读 IRIS 数据服务通讯。IRIS 也已不再维护 rdseed 软件。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import glob
import subprocess

if len(sys.argv) != 2:
    sys.exit("Usage: python {} dirname".format(sys.argv[0]))

dirname = sys.argv[1]

os.chdir(dirname)  # cd 进数据所在目录，以避免处理路径

for seed in glob.glob("*.seed"):
    subprocess.call(['rdseed', '-pdf', seed])

os.chdir("..")

文件合并

SEED 文件的波形数据可能会因为多种原因而出现间断，导致同一个通道会解压出来多个 SAC 文件，因而需要将属于同一个通道的 SAC 数据合并起来。

• 此处使用了新版 merge 命令的语法，要求 SAC 版本大于 v101.6

• merge 命令还有更多选项可以控制数据合并的细节，见命令的语法介绍

• 合并后的数据，以最早的数据段的文件名命名

• 多余的数据段均被删除，以保证每个通道只有一个 SAC 文件

• 由于脚本运行速度比 SAC 运行速度快，因而应先退出 SAC 再删除多余的数据段

下载地址： 2.merge.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import glob
import subprocess

os.putenv("SAC_DISPLAY_COPYRIGHT", '0')

if len(sys.argv) != 2:
    sys.exit("python {} dirname".format(sys.argv[0]))

dirname = sys.argv[1]

os.chdir(dirname)

# 利用 dict 的 key 的不可重复性构建集合:
#     dict 的 key 定义为 NET.STA.LOC.CHN
#     dict 的 value 是文件名与 key 匹配的 SAC 文件数目
sets = {}
for fname in glob.glob("*.SAC"):
    key = '.'.join(fname.split('.')[6:10])
    if key not in sets:
        sets[key] = 1
    else:
        sets[key] += 1

# prepare sac commands
s = "wild echo off \n"
to_del = []
for key, value in sets.items():
    # 仅当 value 大于1时才需要 merge
    if value == 1:
        continue

    print("merge {}: {} traces".format(key, value))
    # Python 的 glob 返回值是乱序的，因而必须 sort
    traces = sorted(glob.glob('.'.join(['*', key, '?', 'SAC'])))

    # 在 SAC 中使用通配符而不是使用 @traces 以避免命令行过长的问题
    # merge 不支持通配符
    s += "r *.{}.?.SAC \n".format(key)    # SAC v101.6 or later
    # 缺数据则补0，重叠则做平均
    s += "merge gap zero overlap average \n"
    s += "w {} \n".format(traces[0])      # 以最早数据段的文件名保存

    to_del.extend(traces[1:])

s += "q \n"
subprocess.Popen(['sac'], stdin=subprocess.PIPE).communicate(s.encode())

# 删除多余的数据段
for file in to_del:
    os.unlink(file)

os.chdir("..")

文件重命名

从 SEED 解压出的 SAC 文件名较长，因而对其重命名以简化。

• SEED 解压出的默认文件名格式为 yyyy.ddd.hh.mm.ss.ffff.NN.SSSSS.LL.CCC.Q.SAC

• 重命名后的文件名为 NET.STA.LOC.CHN.SAC

下载地址： 3.rename.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import glob

os.putenv("SAC_DISPLAY_COPYRIGHT", '0')

if len(sys.argv) != 2:
    sys.exit("Usage: python {} dirname\n".format(sys.argv[0]))

dirname = sys.argv[1]

os.chdir(dirname)
for fname in glob.glob("*.SAC"):
    net, sta, loc, chn = fname.split('.')[6:10]
    os.rename(fname, "{}.{}.{}.{}.SAC".format(net, sta, loc, chn))

os.chdir("..")

添加事件信息

若 SEED 中不包含事件信息，则解压得到的 SAC 文件中也不会包含事件信息。因而需要用户手动添加事件的发震时刻、经纬度、深度和震级信息。

• 输入参数包括：目录名、发震时刻、经度、纬度、深度、震级

• 发震时刻的格式为 yyyy-mm-ddThh:mm:ss.xxx，其中 T 用于分隔日期和时间

下载地址： 4.eventinfo.py

#!/usr/bin/env python
# -*- coding: utf8 -*-

import os
import sys
import datetime
import subprocess

os.putenv("SAC_DISPLAY_COPYRIGHT", '0')

if len(sys.argv) != 7:
    sys.exit("Usage: python {} dirname yyyy-mm-ddThh:mm:ss.xxx "
             "evlo evla evdp mag".format(sys.argv[0]))

dirname, origin, evlo, evla, evdp, mag = sys.argv[1:]

o = datetime.datetime.strptime(origin, '%Y-%m-%dT%H:%M:%S.%f')
# 计算发震日期是一年中的第几天
jday = o.strftime("%j")
# 提取毫秒值
msec = int(o.microsecond / 1000 + 0.5)

os.chdir(dirname)

s = "wild echo off \n"
s += "r *.SAC \n"
s += "synchronize \n"
s += "ch o gmt {} {} {} {} {} {}\n".format(o.year, jday, o.hour,
                                           o.minute, o.second, msec)
s += "ch allt (0 - &1,o&) iztype IO \n"
s += "ch evlo {} evla {} evdp {} mag {} \n".format(evlo, evla, evdp, mag)
s += "wh \n"
s += "q \n"
subprocess.Popen(['sac'], stdin=subprocess.PIPE).communicate(s.encode())

os.chdir("..")

去仪器响应

使用 PZ 文件去仪器响应。若数据的时间跨度太长，在该时间跨度内可能仪器响应会发生变化，因而会存在一个通道有多个 PZ 文件的情况。目前该脚本在遇到这种情况时会自动退出。

下载地址： 5.transfer.py

#!/usr/bin/env python
# -*- coding: utf8 -*-
import os
import sys
import glob
import subprocess

os.putenv("SAC_DISPLAY_COPYRIGHT", "0")

if len(sys.argv) != 2:
    sys.exit("Usage: python {} dirname\n".format(sys.argv[0]))

dirname = sys.argv[1]

os.chdir(dirname)

s = ""
for sacfile in glob.glob("*.SAC"):
    net, sta, loc, chn = sacfile.split('.')[0:4]
    pz = glob.glob("SAC_PZs_{}_{}_{}_{}_*_*".format(net, sta, chn, loc))
    # 暂不考虑多个PZ文件的情况
    if len(pz) != 1:
        sys.exit("PZ file error for {}".format(sacfile))

    s += "r {} \n".format(sacfile)
    s += "rmean; rtr; taper \n"
    s += "trans from pol s {} to none freq 0.01 0.02 5 10\n".format(pz[0])
    s += "mul 1.0e9 \n"
    s += "w over \n"

s += "q \n"
subprocess.Popen(['sac'], stdin=subprocess.PIPE).communicate(s.encode())

os.chdir("..")

分量旋转

将成对的水平分量旋转到大圆路径。

• 检查三分量是否缺失

• 检查 delta 是否相等

• 取三分量中的最大 b 和最小 e 值作为数据窗口，此操作要求三分量的 kzdate 和 kztime 必须相同，这一点在添加事件信息时使用 synchronize 已经实现

• 检查两个水平分量是否正交，若不正交，则无法旋转

下载地址： 6.rotate.py

#!/usr/bin/env python
# -*- coding: utf8 -*-

import os
import sys
import glob
import subprocess

os.putenv("SAC_DISPLAY_COPYRIGHT", '0')

if len(sys.argv) != 2:
    sys.exit("Usage: python %s dirname" % sys.argv[0])

dir = sys.argv[1]
os.chdir(dir)

# 构建 NET.STA.LOC.CH 的集合
sets = set()
for fname in glob.glob("*.SAC"):
    net, sta, loc, chn = fname.split('.')[0:4]
    key = '.'.join([net, sta, loc, chn[0:2]])
    sets.add(key)

p = subprocess.Popen(['sac'], stdin=subprocess.PIPE)
s = ""
for key in sets:
    # 检测Z分量是否存在
    Z = key + "Z.SAC"
    if not os.path.exists(Z):
        print("%s: Vertical component missing!" % key)
        continue

    # 检测水平分量是否存在
    if os.path.exists(key + "E.SAC") and os.path.exists(key + "N.SAC"):
        E = key + "E.SAC"
        N = key + "N.SAC"
    elif os.path.exists(key + "1.SAC") and os.path.exists(key + "2.SAC"):
        E = key + "1.SAC"
        N = key + "2.SAC"
    else:
        print("%s: Horizontal components missings!" % key)
        continue

    # 检查水平分量是否正交
    cmd = 'saclst cmpaz f {}'.format(E).split()
    Ecmpaz = subprocess.check_output(cmd).decode().split()[1]
    cmd = 'saclst cmpaz f {}'.format(N).split()
    Ncmpaz = subprocess.check_output(cmd).decode().split()[1]
    cmpaz_delta = abs(float(Ecmpaz) - float(Ncmpaz))
    if not (abs(cmpaz_delta-90) <= 0.01 or abs(cmpaz_delta-270) <= 0.01):
        print("{}: cmpaz1={}, cmpaz2={} "
              "are not orthogonal!".format(key, Ecmpaz, Ncmpaz))
        continue

    # 检查B, E, DELTA
    cmd = 'saclst b e delta f {}'.format(Z).split()
    Zb, Ze, Zdelta = subprocess.check_output(cmd).decode().split()[1:]
    cmd = 'saclst b e delta f {}'.format(E).split()
    Eb, Ee, Edelta = subprocess.check_output(cmd).decode().split()[1:]
    cmd = 'saclst b e delta f {}'.format(N).split()
    Nb, Ne, Ndelta = subprocess.check_output(cmd).decode().split()[1:]

    if not (float(Zdelta) == float(Edelta) and float(Zdelta) == float(Ndelta)):
        print("%s: delta not equal!" % key)
        continue

    # 获取三分量的最大B和最小E值作为数据窗
    begin = max(float(Zb), float(Eb), float(Nb))
    end = min(float(Ze), float(Ee), float(Ne))

    # 输出文件名为 NET.STA.LOC.[RTZ]
    prefix = key[:-2]
    R, T, Z0 = prefix + '.R', prefix + '.T', prefix + '.Z'

    s += "cut %f %f \n" % (begin, end)
    s += "r %s %s \n" % (E, N)
    s += "rotate to gcp \n"
    s += "w %s %s \n" % (R, T)
    s += "r %s \n" % Z
    s += "w %s \n" % Z0
s += "q \n"
p.communicate(s.encode())

# 删除原文件
for fname in glob.glob("*.SAC"):
    os.unlink(fname)

os.chdir("..")

数据重采样

通常有两种情况下需要对数据进行重采样：

• 原始数据的采样率过高，而实际研究中不需要如此高的采样率，此时，对数据做减采样可以大大减少数据量；

• 原始数据中，不同台站的采样率不同，可能会影响到后期的数据处理，因而需要让所有数据使用统一的采样率；

下面的Python脚本中使用 interpolate 命令将所有数据重采样到相同的采样周期。用户可以在命令行中直接指定要使用的重采样后的采样周期，若命令行中的采样周期指定为0，则以大多数数据所使用的采样周期作为重采样后的采样周期。

下载地址： 7.resample.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import glob
import subprocess

os.putenv("SAC_DISPLAY_COPYRIGHT", "0")

if len(sys.argv) != 3:
    sys.exit("Usage: python {} dirname [delta | 0]".format(sys.argv[0]))
# 若第二个参数为0，则取数据中出现次数最多的周期为采样周期

dirname, delta = sys.argv[1], float(sys.argv[2])

os.chdir(dirname)

if delta == 0:
    # 假定所有数据已旋转到RTZ坐标系，文件名格式为NET.STA.LOC.[RTZ]
    # hash的键为出现了的采样周期，其值为出现的个数
    groups = {}
    for fname in glob.glob("*.*.*.[RTZ]"):
        cmd = 'saclst delta f {}'.format(fname).split()
        delta0 = subprocess.check_output(cmd).decode().split()[1]
        if delta0 not in groups:
            groups[delta0] = 1
        else:
            groups[delta0] += 1
    delta = float(max(groups, key=groups.get))

s = ""
for fname in glob.glob("*.*.*.[RTZ]"):
    cmd = 'saclst delta f {}'.format(fname).split()
    delta0 = float(subprocess.check_output(cmd).decode().split()[1])
    if delta == delta0:  # 不需要重采样
        continue

    s += "r {} \n".format(fname)
    # 用 interpolate 实现减采样或增采样
    # 若是减采样，则需要对数据做低通滤波以防止出现混淆效应
    # 低通滤波时或许需要加上 p 2 以避免滤波引起的相移
    if delta > delta0:
        s += "lp c {} \n".format(0.5/delta)
    s += "interpolate delta {} \n".format(delta)
    s += "w over \n"
s += "q \n"
subprocess.Popen(['sac'], stdin=subprocess.PIPE).communicate(s.encode())

os.chdir("..")