-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathazi_trans.py
125 lines (102 loc) · 4.67 KB
/
azi_trans.py
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
import surfdbase
import copy
# dset = surfdbase.invhdf5('/work1/leon/ALASKA_work/azi_inv_files/azi_20190701_threelay_psimisfit.h5')
dset = surfdbase.invhdf5('/work1/leon/ALASKA_work/azi_inv_files/azi_20190701_lab002.h5')
# dset = surfdbase.invhdf5('/work1/leon/ALASKA_work/azi_inv_files/azi_20190701_trans_0.3_0.75.h5')
# dset = surfdbase.invhdf5('/work1/leon/ALASKA_work/azi_inv_files/azi_20190701_threelay.h5')
#-------------------------
# before inversion
#-------------------------
# dset.read_eik_azi_aniso(inh5fname='/work1/leon/ALASKA_work/hdf5_files/azi_2deg_0.05_20190617.h5')
# -------------------------
# computing kernels
# -------------------------
# dset.compute_kernels_hti(misfit_thresh=5.)
#-------------------------
# LAB
#-------------------------
# dset.read_LAB()
# dset.construct_dvs()
# dset.construct_slab_edge()
# -------------------------
# inversion
# -------------------------
# three layer
# -------------------------
# dset.linear_inv_hti(misfit_thresh=5.)
# dset.construct_hti_model()
# -------------------------
# vpr = dset.compute_kernels_hti(misfit_thresh=5.)
# vpr = dset.compute_kernels_hti(outlon=209., outlat=63.1)
# vpr = dset.linear_inv_hti(outlon=-150.+360., outlat = 65., depth_mid_mantle=80.)
# vpr = dset.linear_inv_hti(outlon=-155.+360., outlat = 63.)
# dset.linear_inv_hti(misfit_thresh=5.)
# dset.construct_hti_model()
# dset.linear_inv_hti_transdimensional(misfit_thresh=5., improve_thresh=.3)
# dset.construct_hti_model_four_lay()
# vpr = dset.linear_inv_hti(outlon=-150.+360., outlat = 65., depth_mid_mantle=80.)
#
# dset.linear_inv_hti_adaptive(misfit_thresh=5., labthresh=80.)
dset.construct_hti_model_three_lay()
# dset.construct_hti_model_four_lay()
#
# # dset.linear_inv_hti(misfit_thresh=5., depth_mid_mantle=80.)
# # dset.construct_hti_model_four_lay()
#
# # d = dset.plot_hti_diff_misfit(inh5fname='/work1/leon/ALASKA_work/azi_inv_files/azi_20190701_fourlay.h5')
# #
cmap = surfdbase.discrete_cmap(6, 'hot_r')
dset.plot_hti(scaled=True, normv=5.,factor=5, gindex=2, datatype='misfit', ampref=.3, plot_data=True, plot_axis=False, cmap=cmap, vmin=0.5, vmax=2.0)
# dset.plot_hti_vel(depth=100., gindex=2, scaled=True, factor=5, ampref=2., normv=1., vmin=4.1, vmax=4.6)
# dset.plot_hti_vel(depth=20., gindex=1, scaled=True, factor=5, ampref=2., normv=1., vmin=3.4, vmax=3.8)
# dset.plot_hti_vel(depth=20., gindex=1, scaled=True, factor=5, ampref=2., normv=2., vmin=3.4, vmax=3.8, ticks=[3.4, 3.5, 3.6, 3.7, 3.8])
# dset.plot_hti_vel(depth=60., gindex=2, scaled=True, factor=5, ampref=2., normv=1., vmin=4.15, vmax=4.55)
# dset.plot_hti_vel(depth=120., gindex=3, scaled=True, factor=5, ampref=2., normv=1., vmin=4.1, vmax=4.6)
##
# LAB
#
# dset.plot_hti(gindex=-1, datatype='labarr', plot_data=True, plot_axis=False)
# dset.plot_paraval(pindex='moho', isthk=True, is_smooth=True, cmap='gist_ncar', vmin=20., vmax=60.0, clabel='Moho Depth (km)')
#
# import matplotlib.pyplot as plt
# import numpy as np
# cmap = plt.cm.RdYlBu
# #
# # # Transparent colours
# # from matplotlib.colors import ListedColormap
# #
# # colA = cmap(np.arange(cmap.N))
# # colA[:,-1] = 0.25 + 0.5 * np.linspace(-1.0, 1.0, cmap.N)**2.0
# #
# # # Create new colormap
# cmap = ListedColormap(colA)
# cmap = surfdbase.discrete_cmap(10, 'RdYlBu')
# cmap = surfdbase.discrete_cmap(10, 'hot_r')
# # # # dset.plot_paraval(pindex='avg_misfit', is_smooth=False, cmap=cmap, vmin=0.0, vmax=2.0, outfname='avg_misfit.txt', clabel='Misfit')
# dset.plot_paraval(pindex='moho', isthk=False, is_smooth=True, cmap=cmap, vmin=25., vmax=45.0, clabel='Crustal thickness (km)')
# dset.plot_paraval(pindex='moho', isthk=True, dtype='std', is_smooth=True, cmap=cmap, vmin=0., vmax=10.0, clabel='Uncertainties of Crustal Thickness (km)')
# # # dset.plot_paraval(pindex='vs_std', is_smooth=False, depth=10., depthavg = 0.)
# #
# dset.plot_crust1( infname='crsthk.xyz', vmin=25., vmax=45., clabel='Crustal thickness (km)', cmap=cmap)
# dset.plot_miller_moho_finer(vmin=25., vmax=45., clabel='Crustal thickness (km)', cmap=cmap)
# dset.convert_to_vts(outdir='outvts', depthavg=-1.)
# # #
# # # import copy
# # # vpr = dset.compute_kernels_hti(misfit_thresh=5., outlon=-155.+360., outlat = 61.3)
# # # vpr2 = copy.deepcopy(vpr)
# # #
# # # vpr.compute_reference_vti(wtype='ray')
# # # vpr2.compute_reference_vti_2(wtype='ray')
# # # v0 = vpr.data.dispR.pvelref.copy()
# # #
# # # vpr.eigkR.bveti[-10:] *= 1.02
# # # vpr2.eigkR.bveti[-10:] *= 1.02
# # # vpr.model.vsv[-10:] *= 1.02
# # # vpr2.model.vsv[-10:] *= 1.02
# # #
# # # v1 = v0 + vpr.eigkR.eti_perturb_vel()
# # # v2 = v0 + vpr2.eigkR.eti_perturb_vel()
# # #
# # # vpr.model.vel2love()
# # # vpr.compute_reference_vti(wtype='ray')
# # # v3 = vpr.data.dispR.pvelref.copy()