Geopy/pre_process.py at master · HouseJaay/Geopy · GitHub

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
from subprocess import Popen,PIPE
from glob import glob
import obspy
import os
from os.path import isfile

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


def trans(filename, respname, writename, f, d):
    p = Popen(['sac'], stdin=PIPE, stdout=PIPE)
    s = ""
    s += "r %s\n" % filename
    s += "decimate %d;decimate %d\n" % (d[0], d[1])
    s += "rmean;rtrend\n"
    s += "transfer from polezero subtype %s to none freq %f %f %f %f\n" % (
        respname, f[0], f[1], f[2], f[3])
    s += "w %s\n" % writename
    s += "q\n"
    r = p.communicate(s.encode())
    #print(r)


def do_trans(sacdir, respdir, writedir):
    saclist = glob(sacdir+'*.BH?')
    print('processing %s' %sacdir)
    for sacpath in saclist:
        sacfile = sacpath.split('/')[-1]
        sta, ch = sacfile.split('.')[1], sacfile.split('.')[-1]
        respname = respdir + "SAC_PZs_X2_%s_%s_00" % (sta, ch)
        writename = writedir + '.'.join(sacfile.split('.')[:-1]) +\
            '.' + ch[-1]
        trans(sacpath, respname, writename, (0.008, 0.01, 3, 4), (5, 2))


def do_rotate(peddir):
    zfiles = glob(peddir+'*.Z')
    for zfile in zfiles:
        common = '.'.join(zfile.split('.')[:-1])
        efile,nfile = common + '.E',common + '.N'
        rfile,tfile = common + '.R',common + '.T'
        if isfile(efile) and isfile(nfile):
            p = Popen(['sac'], stdin=PIPE,stdout=PIPE)
            s = ""
            s += "r %s\n" % efile
            s += "ch cmpinc 90 cmpaz 90;wh\n"
            s += "r %s\n" %nfile
            s += "ch cmpinc 90 cmpaz 0;wh\n"
            s += "r %s %s\n" %(efile, nfile)
            s += "rotate to gcp\n"
            s += "w %s %s\n" % (rfile, tfile)
            s += "q\n"
            r = p.communicate(s.encode())
            #print(r)
        else:
            print('file error %s' %zfile)


if __name__ == '__main__':
    root = '/home/haosj/data/tibet/'
    dirs = os.listdir(root+'ordos/')
    for sta in dirs:
        writedir = root + 'ped/' + sta + '/'
        os.mkdir(writedir)
        do_trans(root+'ordos/'+sta+'/',root+'RESP/RESP/',writedir)
        do_rotate(writedir)