rf2depth_makedata.py 6.11 KB
Newer Older
Mijian Xu's avatar
Mijian Xu committed
1
2
from seispy.rfcorrect import SACStation, psrf2depth, Mod3DPerturbation, psrf_1D_raytracing,\
    psrf_3D_migration, time2depth
Mijian Xu's avatar
Mijian Xu committed
3
4
import numpy as np
from seispy.ccppara import ccppara
Mijian Xu's avatar
Mijian Xu committed
5
from seispy.setuplog import setuplog
Mijian Xu's avatar
Mijian Xu committed
6
from seispy.geo import latlon_from, deg2km, rad2deg
Mijian Xu's avatar
Mijian Xu committed
7
from os.path import join, dirname, exists
Mijian Xu's avatar
Mijian Xu committed
8
from scipy.io import savemat
Mijian Xu's avatar
Mijian Xu committed
9
10
import argparse
import sys
Mijian Xu's avatar
Mijian Xu committed
11
12
13
14


class Station(object):
    def __init__(self, sta_lst):
Mijian Xu's avatar
Mijian Xu committed
15
16
17
        dtype = {'names': ('station', 'evla', 'evlo'), 'formats': ('U20', 'f4', 'f4')}
        self.station, self.stla, self.stlo = np.loadtxt(sta_lst, dtype=dtype, unpack=True, ndmin=1)
        #self.station = [sta.decode() for sta in self.station]
Mijian Xu's avatar
Mijian Xu committed
18
        self.sta_num = self.stla.shape[0]
Mijian Xu's avatar
Mijian Xu committed
19
20
21
22
23
24


def init_mat(sta_num):
    dtype = {'names': ('Station', 'stalat', 'stalon', 'Depthrange', 'events', 'bazi', 'rayp', 'phases', 'moveout_correct',
                      'Piercelat', 'Piercelon', 'StopIndex'),
             'formats': tuple(['O']*12)}
Mijian Xu's avatar
Mijian Xu committed
25
    return np.zeros([1, sta_num], dtype=dtype)
Mijian Xu's avatar
Mijian Xu committed
26
27
28
29
30
31
32
33
34


def _convert_str_mat(instr):
    mat = np.zeros((len(instr), 1), dtype='O')
    for i in range(len(instr)):
        mat[i, 0] = np.array(np.array([instr[i]]), dtype='O')
    return mat


Mijian Xu's avatar
Mijian Xu committed
35
36
37
38
39
40
41
42
43
44
45
def makedata(cpara, velmod3d=None, log=setuplog()):
    if velmod3d is not None:
        if isinstance(velmod3d, str):
            if exists(velmod3d):
                model_3d = np.load(velmod3d)
            else:
                model_3d = None
        else:
            raise ValueError('Path to 3d velocity model should be in str')
    else:
        model_3d = None
Mijian Xu's avatar
Mijian Xu committed
46
    # cpara = ccppara(cfg_file)
Mijian Xu's avatar
Mijian Xu committed
47
    sta_info = Station(cpara.stalist)
Mijian Xu's avatar
Mijian Xu committed
48
    RFdepth = []
Mijian Xu's avatar
Mijian Xu committed
49
    for i in range(sta_info.stla.shape[0]):
Mijian Xu's avatar
Mijian Xu committed
50
        rfdep = {}
Mijian Xu's avatar
Mijian Xu committed
51
        evt_lst = join(cpara.rfpath, sta_info.station[i], sta_info.station[i] + 'finallist.dat')
Mijian Xu's avatar
Mijian Xu committed
52
        stadatar = SACStation(evt_lst, only_r=True)
Mijian Xu's avatar
Mijian Xu committed
53
        log.RF2depthlog.info('the {}th/{} station with {} events'.format(i + 1, sta_info.stla.shape[0], stadatar.ev_num))
Mijian Xu's avatar
Mijian Xu committed
54
55
        piercelat = np.zeros([stadatar.ev_num, cpara.depth_axis.shape[0]])
        piercelon = np.zeros([stadatar.ev_num, cpara.depth_axis.shape[0]])
Mijian Xu's avatar
Mijian Xu committed
56
57
        PS_RFdepth, end_index, x_s, x_p = psrf2depth(stadatar, cpara.depth_axis, stadatar.sampling, stadatar.shift, cpara.velmod, 
                                             velmod_3d=model_3d, srayp=cpara.rayp_lib)
Mijian Xu's avatar
Mijian Xu committed
58
59
        for j in range(stadatar.ev_num):
            piercelat[j], piercelon[j] = latlon_from(sta_info.stla[i], sta_info.stlo[i],
Mijian Xu's avatar
Mijian Xu committed
60
                                                     stadatar.bazi[j], rad2deg(x_s[j]))
Mijian Xu's avatar
Mijian Xu committed
61
62
63
        rfdep['Station'] = sta_info.station[i]
        rfdep['stalat'] = sta_info.stla[i]
        rfdep['stalon'] = sta_info.stlo[i]
Mijian Xu's avatar
Mijian Xu committed
64
65
        # rfdep['Depthrange'] = cpara.depth_axis
        # rfdep['events'] = _convert_str_mat(stadatar.event)
Mijian Xu's avatar
Mijian Xu committed
66
67
        rfdep['bazi'] = stadatar.bazi
        rfdep['rayp'] = stadatar.rayp
Mijian Xu's avatar
Mijian Xu committed
68
        # rfdep['phases'] = _convert_str_mat(stadatar.phase)
Mijian Xu's avatar
Mijian Xu committed
69
70
71
        rfdep['moveout_correct'] = PS_RFdepth
        rfdep['Piercelat'] = piercelat
        rfdep['Piercelon'] = piercelon
Mijian Xu's avatar
Mijian Xu committed
72
73
        rfdep['StopIndex'] = end_index
        RFdepth.append(rfdep)
Mijian Xu's avatar
Mijian Xu committed
74
75
    savemat(cpara.depthdat, {'RFdepth': RFdepth})
    # np.save(cpara.depthdat, RFdepth)
Mijian Xu's avatar
Mijian Xu committed
76
77


Mijian Xu's avatar
Mijian Xu committed
78
79
80
def makedata3d(cpara, velmod3d, log=setuplog()):
    mod3d = Mod3DPerturbation(velmod3d, cpara.depth_axis, velmod=cpara.velmod)
    sta_info = Station(cpara.stalist)
Mijian Xu's avatar
Mijian Xu committed
81
82
83
84
    if cpara.rayp_lib is not None:
        srayp = np.load(cpara.rayp_lib)
    else:
        srayp = None
Mijian Xu's avatar
Mijian Xu committed
85
86
87
88
89
90
    RFdepth = []
    for i in range(sta_info.stla.shape[0]):
        rfdep = {}
        evt_lst = join(cpara.rfpath, sta_info.station[i], sta_info.station[i] + 'finallist.dat')
        stadatar = SACStation(evt_lst, only_r=True)
        log.RF2depthlog.info('the {}th/{} station with {} events'.format(i + 1, sta_info.stla.shape[0], stadatar.ev_num))
91
        pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, tpds = psrf_1D_raytracing(stadatar,
Mijian Xu's avatar
Mijian Xu committed
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
                                                                                                cpara.depth_axis, srayp=srayp)
        newtpds = psrf_3D_migration(pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p,
                                    tpds, cpara.depth_axis, mod3d)
        amp3d, end_index = time2depth(stadatar, cpara.depth_axis, newtpds)
        rfdep['Station'] = sta_info.station[i]
        rfdep['stalat'] = sta_info.stla[i]
        rfdep['stalon'] = sta_info.stlo[i]
        # rfdep['Depthrange'] = cpara.depth_axis
        # rfdep['events'] = _convert_str_mat(stadatar.event)
        rfdep['bazi'] = stadatar.bazi
        rfdep['rayp'] = stadatar.rayp
        # rfdep['phases'] = _convert_str_mat(stadatar.phase)
        rfdep['moveout_correct'] = amp3d
        rfdep['Piercelat'] = pplat_s
        rfdep['Piercelon'] = pplon_s
        rfdep['StopIndex'] = end_index
        RFdepth.append(rfdep)
Mijian Xu's avatar
Mijian Xu committed
109
110
111
112
113
114
    try:
        savemat(cpara.depthdat, {'RFdepth': RFdepth})
    except FileNotFoundError:
        log.RF2depthlog.warning('No such file or directory: {}'.format(dirname(cpara.depthdat)))
        rfdep_path = input('Enter a exist path:')
        savemat(rfdep_path, {'RFdepth': RFdepth})
Mijian Xu's avatar
Mijian Xu committed
115
116


Mijian Xu's avatar
Mijian Xu committed
117
118
def rf2depth():
    parser = argparse.ArgumentParser(description="Convert Ps RF to depth axis")
Mijian Xu's avatar
Mijian Xu committed
119
120
    parser.add_argument('-d', help='Path to 3d vel model in npz file for moveout correcting', type=str, default='')
    parser.add_argument('-r', help='Path to 3d vel model in npz file for 1D ray tracing', type=str, default='')
Mijian Xu's avatar
Mijian Xu committed
121
122
123
124
125
    parser.add_argument('cfg_file', type=str, help='Path to configure file')
    arg = parser.parse_args()
    if len(sys.argv) == 1:
        parser.print_help()
        sys.exit(1)
Mijian Xu's avatar
Mijian Xu committed
126
    cpara = ccppara(arg.cfg_file)
Mijian Xu's avatar
Mijian Xu committed
127
128
129
130
131
132
    if arg.d != '' and arg.r != '':
        raise ValueError('Specify only 1 argument in \'-d\' and \'-r\'')
    elif arg.d != '' and arg.r == '':
        makedata3d(cpara, arg.d)
    elif arg.d == '' and arg.r != '':
        makedata(cpara, arg.r)
Mijian Xu's avatar
Mijian Xu committed
133
    else:
Mijian Xu's avatar
Mijian Xu committed
134
        makedata(cpara)
Mijian Xu's avatar
Mijian Xu committed
135
136


Mijian Xu's avatar
Mijian Xu committed
137
if __name__ == '__main__':
Mijian Xu's avatar
Mijian Xu committed
138
    cfg_file = '/Users/xumj/Researches/Tibet_MTZ/process/paraCCP.cfg'
Mijian Xu's avatar
Mijian Xu committed
139
    vel3d_file = '/Users/xumj/Researches/Tibet_MTZ/models/GYPSUM.npz'
Mijian Xu's avatar
Mijian Xu committed
140
    cpara = ccppara(cfg_file)
Mijian Xu's avatar
Mijian Xu committed
141
    makedata3d(cpara, vel3d_file)