Huston PlotΒΆ
The example_path need to de changed to your path, and run this notebook.
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3import os
4import numpy as np
5import matplotlib.pyplot as plt
6from matplotlib import cm
7from MCMTpy import MomentTensor as MTpy
8
9# Helper function
10#----------------------------------------------------#
11#%% 4. plot_Hudson_points
12def plot_Hudson_points(FM):
13 fig1, ax1 = plt.subplots(1 ,1, dpi=800)
14 plt.rc('font',family='Times New Roman')
15 MTpy.Hudson_plot(ax=ax1)
16 for i in range(0,len(FM)):
17 MT = MTpy.MTensor(FM[i,:])
18 M=MT.mt
19
20 k,T = MTpy.M2kT_space(M)
21 U,V = MTpy.kT2UV_space(k,T)
22
23 map_vir = cm.get_cmap(name='YlGn')
24 colors = map_vir(i/len(FM))
25
26 ax1.scatter(U,V, color=colors,marker='o', s=5)
27
28
29 position=fig1.add_axes([0.85, 0.15, 0.01, 0.5])
30 font_colorbar = {'family' : 'Times New Roman',
31 'color' : 'black',
32 'weight' : 'normal',
33 'size' : 6,
34 }
35 sm = cm.ScalarMappable(cmap=map_vir)
36 sm.set_array(np.arange(0,len(FM)+1))
37 cb=plt.colorbar(sm,cax=position,orientation='vertical')
38 cb.set_label('Sample',fontdict=font_colorbar)
39
40 ax1.set_xlim(-4/3-0.1, 4/3+0.3)
41 ax1.set_ylim(-1-0.1, 1+0.1)
42
43 return fig1
44
45
46#%%------------------------------------------------------------#
47## 0. set path
48example_path = '/Users/yf/3.Project/8.MCMTpy/MCMTpy-master/data/example_yunnan'
49
50
51#------------------------------------------------------------#
52# we expect no parameters need to be changed below
53#------------------------------------------------------------#
54FM_path=os.path.join(example_path,"YN.202105212148_Inv/mt_inv/Output_YN.202105212148_mt/rank_0_output/chain_0_FM_accept_all")
55allfiles_path = os.path.join(example_path,'YN.202105212148_Inv/YN.202105212148_raw/*.SAC')
56## 1. read FM
57N=6
58FM_all = np.loadtxt(FM_path)
59FM_raw = FM_all[0:,1:N+1]
60Mxx_np = FM_raw[:,0]
61Mxy_np = FM_raw[:,1]
62Mxz_np = FM_raw[:,2]
63Myy_np = FM_raw[:,3]
64Myz_np = FM_raw[:,4]
65Mzz_np = FM_raw[:,5]
66FM = np.vstack((Mxx_np, Myy_np, Mzz_np, Mxy_np, Mxz_np, Myz_np)).T
67
68## 2. plot Hudson
69fig = plot_Hudson_points(FM)
70figurename=os.path.join('./S2_figure/Hudson.pdf')
71plt.savefig(figurename,dpi=800, format="pdf")
