-
Notifications
You must be signed in to change notification settings - Fork 16
How to use SeisTomoPy class
Go in the downloaded directory and launch iPython. The first thing to do is to import some packages, including SeisTomoPy, so that you will be able to use the various SeisTomoPy functions.
import SeisTomoPy
import matplotlib.pyplot as plt
import numpy as np
The user can generate cross sections through Vp, Vs or ρ variations anywhere on Earth in some default tomographic models provided with the tool using SeisTomoPy.cross_section_plot. This generates directly the plot of the cross section. However, the user could also be interested in getting the values of the cross section and then perform some further calculations. This can also be achieved running now SeisTomoPy.cross_section. The latter function returns as outputs the matrix containing the cross section (Z) as well as the angle (th) and radius vectors (r) that the user would need if he desires to perform further calculations with this cross section.
Below is an example of how to get a cross section beneath Africa using both class. Pink circles show hotspot locations, white circles denote earthquake locations and green, red and yellow circles the starting, ending and mid-point, respectively, along the profile.
#Setting parameters
#Model to be plotted
#Choose between:
#SEISGLOB2, S40RTS, SEMUCBWM1, S362WMANIM, SEISGLOB1, SP12RTS, SGLOBE, 3D2016, MYMODEL
model = 'SEISGLOB2'
#Parameter to be plotted
#Choose between: VS, VP, RHO
para = 'VS'
#Latitude of the starting point of the cross section
elat = -60
#Longitude of the starting point of the cross section
elon = -49
#Latitude of the ending point of the cross section
slat = 60
#Longitude of the ending point of the cross section
slon = 119
#Depth of the cross section
depth = 2890
#Spherical harmonic degrees to be used
NSmax = 40
#Maximal velocity perturbations for the colorbar
Vmax = 2
#Running cross_section_plot
SeisTomoPy.cross_section_plot(model,para,elat,elon,slat,slon,depth,NSmax,Vmax)
#Running cross_section
Z, th, r = SeisTomoPy.cross_section(model,para,elat,elon,slat,slon,depth,NSmax)
The user can create maps at a given depth for the whole globe using SeisTomoPy.tomomap_plot. This generates directly the plot of the map. However, the user could also be interested in getting the values of the map and then perform some further calculations. This can be achieved using SeisTomoPy.tomomap.
Below is an example of how to obtain a map for the tomographic model SEISGLOB2 at 50 km depth using both class. Pink circles show hotspot locations and grey lines the plate boundaries.
# Setting parameters
# Model to be plotted
model = '3D2016'
# Parameter to be plotted
para = 'VS'
# Depth of the map to be plotted
depth = 50
# Spherical harmonic degrees to be used
NSmax = 60
# Central longitude
# default value 140 degrees
lon0 = 140
# Maximal velocity perturbations for the color bar
Vmax = 2
# Running tomomap_plot
SeisTomoPy.tomomap_plot(model,para,depth,NSmax,Vmax,lon0)
# Changing the central longitude
lon0 = 0
SeisTomoPy.tomomap_plot(model,para,depth,NSmax,Vmax,lon0)
# Running tomomap
Z, lat, lon = SeisTomoPy.tomomap(model,para,depth,NSmax)
You can also change NSmax to a value below 60, you will then obtain a filtered image of the model.
# Filtering the model
NSmax = 18
SeisTomoPy.tomomap_plot(model,para,depth,NSmax,Vmax)
The user can compute the amplitude spectrum, S(X,z,l), for a given model at a given depth, z, for a given parameter, X, and for a certain spherical harmonic degree, l, using SeisTomoPy.spectrum and SeisTomoPy.spectrum_fromfile. The first one compute the spectrum for the models included by default in SeisTomoPy while the second one enables the user to compute the spectrum in any model he would like as long as he can provide the required input file. Please refer to section input files to have details about these files.
Below is an example of obtained spectrum for various models at 520 km depth, for parameter Vs and up to spherical harmonic degree 40 and then how to get the spectrum for any model not included in SeisTomoPy.
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print(' Example of spectrum computed in the default tomographic models ')
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
# Setting parameters
# Models to be used for computing the spectrum
model1 = 'SEISGLOB2'
model2 = 'SGLOBE'
model3 = 'S40RTS'
model4 = 'SEMUCBWM1'
model5 = 'S362WMANIM'
# Depth at which the spectrum is computed
depth = 520
# Parameter used for computing the spectrum
para = 'VS'
# Maximum spherical harmonic degree up to which the spectrum is computed
NSmax = 40
# Running spectrum
sp1 = SeisTomoPy.spectrum(model1,para,depth,NSmax)
sp2 = SeisTomoPy.spectrum(model2,para,depth,NSmax)
sp3 = SeisTomoPy.spectrum(model3,para,depth,NSmax)
sp4 = SeisTomoPy.spectrum(model4,para,depth,NSmax)
sp5 = SeisTomoPy.spectrum(model5,para,depth,NSmax)
# Plotting the results
plt.plot(sp1[:,0], sp1[:,1]/np.amax(sp1[:,1]), linewidth=1.0,color="red",marker="d",markeredgecolor="k", label=model1)
plt.plot(sp2[:,0], sp2[:,1]/np.amax(sp2[:,1]), linewidth=1.0,color="k", marker="d",markeredgecolor="k", label=model2)
plt.plot(sp3[:,0], sp3[:,1]/np.amax(sp3[:,1]), linewidth=1.0,color="blue",marker="d",markeredgecolor="k", label=model3)
plt.plot(sp4[:,0], sp4[:,1]/np.amax(sp4[:,1]), linewidth=1.0,color="dodgerblue",marker="d",markeredgecolor="k", label=model4)
plt.plot(sp5[:,0], sp5[:,1]/np.amax(sp5[:,1]), linewidth=1.0,color="cyan",marker="d",markeredgecolor="k", label=model5)
plt.legend(bbox_to_anchor=(1.45, 1))
plt.xlabel("Harmonic degree l")
plt.ylabel("Spectrum amplitude")
plt.xlim([0.5, NSmax+0.5])
plt.grid(color='k', linestyle='--', linewidth=0.5)
plt.rcParams.update({'font.size': 12})
plt.show()
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print(' Example of spectrum computed in another model')
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
# Running spectrum_fromfile
filename = 'SeisTomoPy_notebook/files/input_file_spectrum.xyz'
sp_fromfile = SeisTomoPy.spectrum_fromfile(filename,NSmax)
# Plotting the results
plt.plot(sp_fromfile[:,0], sp_fromfile[:,1]/np.amax(sp_fromfile[:,1]),linewidth=1.0,color="orange",marker="d",markeredgecolor="k", label='Spectrum from a file')
plt.legend(bbox_to_anchor=(1.55, 1))
plt.xlabel("Harmonic degree l")
plt.ylabel("Spectrum amplitude")
plt.xlim([0.5, NSmax+0.5])
plt.grid(color='k', linestyle='--',linewidth=0.5)
plt.rcParams.update({'font.size': 12})
plt.show()
The user can compute the correlation between two tomographic models, (1) and (2), using SeisTomoPy.correlation, SeisTomoPy.correlation_fromfile and SeisTomoPy.correlation_fromfile2. The first one compute the correlation between two models chosen from the included ones by default in SeisTomoPy, while the second and third ones enable the user to compute the correlation between either any model with one of the default onee or between any other models that the user would have. In order to use the two last functions the user thus must provide the required input file. Please refer to section input files to have details about these files.
Below is an example of the correlation between SEISGLOB2 and S40RTS computed at 520 km depth, for parameter Vs and up to spherical harmonic degree 40 and then how to get the correlation for any model not included in SeisTomoPy.
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print('Example of correlation computed between the default tomographic models')
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
# Setting parameters
# Models to be used for computing the correlation
model1 = 'SEISGLOB2'
model2 = 'S40RTS'
# Depth at which the correlation is computed
depth = 520
# Parameters used for computing the correlation
para1 = 'VS'
para2 = 'VS'
# Maximum spherical harmonic degree up to which the correlation is computed
NSmax = 40
# Running correlation
corr = SeisTomoPy.correlation(model1,model2,depth,depth,para1,para2,NSmax)
plt.plot(corr[:,0],corr[:,1],linewidth=1.0,color="grey",marker="d", label=model1 + '-' + model2)
conf66 = np.loadtxt('conf66.dat')
conf90 = np.loadtxt('conf90.dat')
conf95 = np.loadtxt('conf95.dat')
plt.plot(conf66[:,0], conf66[:,1], linewidth=1.0,color="silver", ls='--', label='Confidence level 66%')
plt.plot(conf90[:,0], conf90[:,1], linewidth=1.0,color="gray", ls='--', label='Confidence level 90%')
plt.plot(conf95[:,0], conf95[:,1], linewidth=1.0,color="black", ls='--', label='Confidence level 95%')
plt.plot(conf66[:,0], -conf66[:,1], linewidth=1.0,color="silver", ls='--')
plt.plot(conf90[:,0], -conf90[:,1], linewidth=1.0,color="gray", ls='--')
plt.plot(conf95[:,0], -conf95[:,1], linewidth=1.0,color="black", ls='--')
plt.xlabel("Harmonic degree l")
plt.ylabel("Correlation")
plt.xlim([0.5, NSmax+0.5])
plt.ylim([-1, 1])
plt.legend(bbox_to_anchor=(1.6, 1))
plt.grid(color='k', linestyle='--', linewidth=0.5)
plt.rcParams.update({'font.size': 12})
plt.show()
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print(' Example of correlation using other models ')
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
# Setting new parameters
# Models to be used for computing the correlation
model1 = 'S40RTS'
# Depth at which the correlation is computed
depth1 = 2800
# Parameters used for computing the correlation
para1 = 'VS'
# Maximum spherical harmonic degree up to which the correlation is computed
NSmax = 40
# Running correlation_fromfile
filename = 'SeisTomoPy_notebook/files/input_file_spectrum.xyz'
corr_fromfile = SeisTomoPy.correlation_fromfile(model1,depth1,para1,filename,NSmax)
# Plotting the results
plt.plot(corr_fromfile[:,0], corr_fromfile[:,1]/np.amax(corr_fromfile[:,1]), linewidth=1.0,color="grey",marker="d",label=model1 + '-from file')
plt.legend(bbox_to_anchor=(1.55, 1))
plt.xlabel("Harmonic degree l")
plt.ylabel("Spectrum amplitude")
plt.xlim([0.5, NSmax+0.5])
plt.grid(color='k', linestyle='--', linewidth=0.5)
plt.rcParams.update({'font.size': 12})
plt.show()
# Running correlation_fromfile2
filename1 = 'SeisTomoPy_notebook/files/input_file_spectrum.xyz'
filename2 = 'SeisTomoPy_notebook/files/input_file_spectrum2.xyz'
corr_fromfile2 = SeisTomoPy.correlation_fromfile2(filename1,filename2,NSmax)
# Plotting the results
plt.plot(corr_fromfile2[:,0], corr_fromfile2[:,1]/np.amax(corr_fromfile2[:,1]), linewidth=1.0,color="grey",marker="d",label='from file1-from file2')
plt.legend(bbox_to_anchor=(1.55, 1))
plt.xlabel("Harmonic degree l")
plt.ylabel("Spectrum amplitude")
plt.xlim([0.5, NSmax+0.5])
plt.grid(color='k', linestyle='--', linewidth=0.5)
plt.rcParams.update({'font.size': 12})
plt.show()
The user may want to check which seismic structures of the mantle are sampled by seismic waves. To do so he can use SeisTomoPy.path_plot and SeisTomoPy.path_plot_fromfile to display seismic wave paths on top of cross sections made in the desired tomographic model coming from either one of the default ones included in SeisTomoPy or that the user provides. Below is an example with tomographic model SEISGLOB2, parameter Vs and seismic phases S, ScS, PKP, PKiKP and Sdiff. If it happens that the seismic phase does not exist for the distance range between the earthquake and the station, then it will simply be ignored.
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print('Two examples of path plots using the default tomographic models')
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
# Setting parameters
# Model to be plotted
model = 'SEISGLOB2'
# Parameter to be plotted
para = 'VS'
# Latitude of the starting point of the cross section
elat = -60
# Longitude of the starting point of the cross section
elon = -49
# Latitude of the ending point of the cross section
slat = 60
# Longitude of the ending point of the cross section
slon = 119
# Maximal velocity perturbations for the color bar
Vmax = 2
# List of seismic phases
phlist = 'S'
# Position of stations and event
EVT = np.loadtxt('SeisTomoPy_notebook/files/event.xy')
STA = np.loadtxt('SeisTomoPy_notebook/files/station.xy')
# Running path_plot
SeisTomoPy.path_plot(model,para,Vmax,elat,elon,slat,slon,EVT,STA,phlist)
# List of seismic phases
phlist = 'S ScS PKP PKiKP Sdiff'
# Running path_plot
SeisTomoPy.path_plot(model,para,Vmax,elat,elon,slat,slon,EVT,STA,phlist)
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
print(' An example of path plot using a model file ')
print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
# Setting new parameters
# Latitude of the starting point of the cross section
elat = -60
# Longitude of the starting point of the cross section
elon = -49
# Latitude of the ending point of the cross section
slat = 60
# Longitude of the ending point of the cross section
slon = 130
# Parameter to be plotted
para = 'VS'
# Angular size of the cross section
width = 180
# Sampling rate in theta direction
width_step = 1
# Theta vector
th = np.arange(90-(width/2), 90+(width/2)+width_step, width_step)
# Starting and ending radius of the cross section
r_min = 3481
r_max = 6281
# Sampling rate in radius direction
r_step = 50
# Radius vector
r = np.arange(r_min, r_max+r_step, r_step)
# List of seismic phases
phlist = 'Pdiff'
# Running path_plot_from_file
filename = 'SeisTomoPy_notebook/files/input_file_path.xyz'
SeisTomoPy.path_plot_fromfile(filename,th,r,para,elat,elon,slat,slon,Vmax,EVT,STA,phlist)
Finally, the user can compute the travel time delays with respect to a reference model through a given tomographic model, dt2D_TomoPy, for any given seismic phases and for any combinations of source and receiver provided by the user using SeisTomoPy.get_travel_time. We provide an example of delays computed in model SEISGLOB2 for S and ScS seismic phases in the same region.
# Setting parameters
# Model to be plotted
model = 'SEISGLOB2'
# Latitude of the event
elat = -21
# Longitude of the event
elon = -179
# Depth of the event
edepth = 610
# Position of stations
STA = np.loadtxt('SeisTomoPy_notebook/files/lat_lon_ttstation.txt')
tt2D = np.zeros(len(STA))
dtt2D = np.zeros(len(STA))
ttREF = np.zeros(len(STA))
dist = np.arange(40,74,1)
degmin = 40
degmax = 73
# Running TimePy routinely
for k in range(len(STA)):
# List of seismic phases
List = ['S']
tt2D[k], dtt2D[k], ttREF[k], phase_name = SeisTomoPy.get_travel_time(model,elat,elon,edepth,STA[k,0],STA[k,1],List)
file_str = str(k) + ' ' + 'S ' + str(tt2D[k]) + ' ' + str(tt2D[k]+ttREF[k]) + ' ' + str(dtt2D[k]) + '\n'
print(file_str)
plt.plot(dist,tt2D,marker="d",linewidth=1.0,color="blue", ls='--',label = 'S')
for k in range(len(STA)):
# List of seismic phases
List = ['ScS']
tt2D[k], dtt2D[k], ttREF[k], phase_name = SeisTomoPy.get_travel_time(model,elat,elon,edepth,STA[k,0],STA[k,1],List)
file_str = str(k) + ' ' + 'ScS ' + str(tt2D[k]) + ' ' + str(tt2D[k]+ttREF[k]) + ' ' + str(dtt2D[k]) + '\n'
print(file_str)
# Plotting the results
plt.plot(dist,tt2D,marker="d",linewidth=1.0,color="red", ls='-',label = 'ScS')
plt.xlabel("Epicentral Distance (degrees)")
plt.ylabel("dt (s)")
plt.xlim([degmin-0.5, degmax+0.5])
plt.grid(color='k', linestyle='--', linewidth=0.5)
plt.rcParams.update({'font.size': 12})
plt.legend(bbox_to_anchor=(1.6, 1))
plt.show()