In [1]:
import pandas as pd
import os
import numpy as np
import matplotlib.pylab as plt
from scipy import optimize
from scipy import stats
from scipy import interpolate
from matplotlib import cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.stats import pearsonr
# stats.binned_statistic() method
from scipy import stats
import warnings
warnings.filterwarnings("ignore")
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
%matplotlib inline
def fullprint(df):
with pd.option_context('display.max_rows', None, 'display.max_columns', None): # more options can be specified also
print(df)
file = os.path.join(os.getcwd(),'rundir4983_stitched_rakes_dips_lengths_m6.5.csv')
rsqsim = pd.read_csv(file,sep=',',skipinitialspace=True,low_memory=False,skiprows=None)
rsqsim.dropna(inplace=True)
rsqsim=rsqsim.rename(columns={"Event ID":"Eventid",
"Occurrence Time (s)": "OccurT",
"Magnitude": "M",
"Moment (N-m)":"Moment",
"Area (m^2)": "A",
"Average Slip (m)":"AD",
"Average Element Slip Rate (m/yr)":"SR",
"Hypocenter Latitude":"Hlat",
"Hypocenter Longitude":"Hlon",
"Centroid Latitude":"Clat",
"Centroid Longitude":"Clon",
"Hypocenter Depth (km)":"HD",
"Centroid Depth (km)":"CD",
"Upper Depth (km)":"TD",
"Lower Depth (km)":"BD",
"Average Rake":"ARake",
"Average Absolute Rake Deviation":"AdRake",
"Max Absolute Rake Deviation":"MdRake",
"Average Dip":"ADip",
"Average Absolute Dip Deviation":"AdDip",
"Max Absolute Dip Deviation":"MdDip",
"Full Mapped Subsection Length (km)":"FML",
"Full Slipped Length (km)":"FSL",
"Mid-Seismogenic Slipped Length (km)":"MSL",
"Surface Slipped Length (km)":"SSL"})
rsqsim['A']=rsqsim['A']/1e6
rsqsim['L']=rsqsim['FSL']
rsqsim['W']=rsqsim['A']/rsqsim['L']
rsqsim['Aratio']=rsqsim['L']/rsqsim['W']
rsqsim['SF']=rsqsim['SSL']/rsqsim['L']
rsqsim['SR']=rsqsim['SR']*1e3 #convert to mm/yr
rsqsim['TW']=rsqsim['TD']/np.sin(rsqsim['ADip']*np.pi/180.)
rsqsim['BW']=rsqsim['BD']/np.sin(rsqsim['ADip']*np.pi/180.)
def chinnery_strdrop(LE,AD,WE=15):
LH=LE/2
a=np.sqrt(LH**2+WE**2)
C1=2*LH/(a*WE)+3/LH-(LH*(3*a+4*WE))/(a*(a+WE)**2)
G=30*1e3
return G*AD*1e-3*C1/2/np.pi
rsqsim['SDP']=chinnery_strdrop(rsqsim['L'].to_numpy(),rsqsim['AD'].to_numpy(),rsqsim['W'].to_numpy())
rsqsim['FM']='Oblique'
rsqsim.loc[(rsqsim['ARake'] >= -180) & (rsqsim['ARake'] <= -170), 'FM'] = 'RSS'
rsqsim.loc[(rsqsim['ARake'] <= 180) & (rsqsim['ARake'] >= 170), 'FM'] = 'RSS'
rsqsim.loc[(rsqsim['ARake'] <= 10) & (rsqsim['ARake'] >= -10), 'FM'] = 'LSS'
rsqsim.loc[(rsqsim['ARake'] <= 100) & (rsqsim['ARake'] >= 80), 'FM'] = 'RV'
rsqsim.loc[(rsqsim['ARake'] <= -80) & (rsqsim['ARake'] >= -100), 'FM'] = 'NM'
In [2]:
import matplotlib.colors as mcolors
def make_colormap(seq,N=256):
"""Return a LinearSegmentedColormap
seq: a sequence of floats and RGB-tuples. The floats should be increasing
and in the interval (0,1).
"""
seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3]
cdict = {'red': [], 'green': [], 'blue': []}
for i, item in enumerate(seq):
if isinstance(item, float):
r1, g1, b1 = seq[i - 1]
r2, g2, b2 = seq[i + 1]
cdict['red'].append([item, r1, r2])
cdict['green'].append([item, g1, g2])
cdict['blue'].append([item, b1, b2])
return mcolors.LinearSegmentedColormap('CustomMap', cdict, N=N)
c = mcolors.ColorConverter().to_rgb
rvb = make_colormap(
[c('yellow'), c('brown'),0.25,c('brown'),c('gold'), 0.5, c('gold'),c('violet'), 0.75,c('violet'), c('dodgerblue')],10)
For strike-slip events¶
In [3]:
fig = plt.figure( figsize=(15,5))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='RSS') | (rsqsim['FM']=='LSS')
axe[0].hist(rsqsim.loc[condition,'ADip'],bins=40,density=True);
axe[1].hist(rsqsim.loc[condition,'AdDip'],bins=40,density=True);
axe[2].hist(rsqsim.loc[condition,'MdDip'],bins=40,density=True);
axe[3].hist(rsqsim.loc[condition,'ARake'],bins=40,density=True);
axe[4].hist(rsqsim.loc[condition,'AdRake'],bins=40,density=True);
axe[5].hist(rsqsim.loc[condition,'MdRake'],bins=40,density=True);
In [4]:
cmap=plt.get_cmap('RdPu',10)
sr_range=[0.01,1,10,40]
fig = plt.figure( figsize=(15,12))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=((rsqsim['FM']=='RSS') | (rsqsim['FM']=='LSS')) & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
x=rsqsim.loc[condition,'SR']
binx = np.logspace(-2,2,50)
field=['M','SDP','L','W','A','Aratio']
binyarray=[np.linspace(6.5,8.1,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[False,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('SR (mm/yr)')
plt.tight_layout()
In [5]:
cmap=plt.get_cmap('RdPu',10)
sr_range=[0.01,1,10,40]
fig = plt.figure( figsize=(15,12))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=((rsqsim['FM']=='RSS') | (rsqsim['FM']=='LSS')) & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
x=rsqsim.loc[condition,'M']
binx = np.linspace(6.5,8.1,50)
field=['SR','SDP','L','W','A','Aratio']
binyarray=[np.logspace(-2,2,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[True,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
# axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('Magnitude')
plt.tight_layout()
In [6]:
cmap=plt.get_cmap('RdPu',10)
c=['darkblue','green','orange']
sr_range=[0.01,1,10,40]
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
fig = plt.figure( figsize=(15,20))
axe=fig.subplots(6, 4)
condition=((rsqsim['FM']=='RSS') | (rsqsim['FM']=='LSS')) & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
binx = np.linspace(6.5,8.3,50)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[True,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
xall=rsqsim.loc[condition,'M']
yall=rsqsim.loc[condition,field[iff]]
if logy[iff]:
yall = np.log10(yall)
meanall,bin_edges,_=stats.binned_statistic(xall, yall, statistic='mean', bins=nbin)
if logy[iff]:
meanall=10**meanall
for i in range(len(sr_range)-1):
condition2=condition & ((rsqsim['SR']>=sr_range[i]) & (rsqsim['SR']<sr_range[i+1]))
y=rsqsim.loc[condition2,field[iff]]
x=rsqsim.loc[condition2,'M']
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff,i].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
# ax1_divider = make_axes_locatable(axe[iff,i])
# cax1 = ax1_divider.append_axes("right", size="4%", pad="-10%")
# cb1 = fig.colorbar(im1, cax=cax1,orientation='vertical')
# cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff,i].set_yscale('log')
axe[iff,i].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,i].grid(which='minor', color='silver', linestyle='--')
axe[iff,i].set_axisbelow(True)
axe[iff,i].set_ylabel(field[iff],fontsize='large')
axe[iff,i].set_xlabel('Magnitude')
x = x.to_numpy()
y = y.to_numpy()
if logy[iff]:
y = np.log10(y)
if(len(y)>0):
mean,bin_edges,_=stats.binned_statistic(x, y, statistic='mean', bins=nbin)
mean_centers = (bin_edges[1:] + bin_edges[0:-1])/2
if logy[iff]:
mean=10**mean
if plotratio:
axe[iff,3].plot(mean_centers,mean/meanall,'-',color=c[i],linewidth=4,zorder=10,label='Mean for {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]) )
axe[iff,3].set_ylim([0.1,10])
axe[iff,3].set_yscale('log')
else:
if(len(y)>0):
axe[iff,3].plot(mean_centers,mean,'-',color=c[i],linewidth=4,zorder=10,label='Mean for {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]) )
if logy[iff]:
axe[iff,3].set_yscale('log')
axe[iff,3].set_ylim([biny.min(),biny.max()])
if iff==0:
axe[iff,3].legend(loc='best',fontsize='medium')
axe[iff,3].set_xlim([binx.min(),binx.max()])
axe[iff,3].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,3].grid(which='minor', color='silver', linestyle='--')
axe[iff,3].set_axisbelow(True)
plt.tight_layout()
In [7]:
cmap=plt.get_cmap('RdPu',10)
fig = plt.figure( figsize=(10,10))
axe=[]
axe.append(plt.subplot2grid((1,1), (0,0), rowspan=1, colspan=1))
condition=((rsqsim['FM']=='RSS') | (rsqsim['FM']=='LSS')) & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
x=rsqsim.loc[condition,'SR']
binx = np.logspace(-2,2,50)
field=['SDP',]
binyarray=[
np.logspace(0,1,50),
]
logy=[True,]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap,alpha=0.8,rasterized=True)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('SR (mm/yr)')
nbin=5
rawx=x.to_numpy()
rawy=y.to_numpy()
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(np.log(rsqsim.loc[condition,'SR'].to_numpy()).reshape((-1, 1)), np.log(rsqsim.loc[condition,'SDP'].to_numpy()).reshape((-1, 1)))
a=model.intercept_
b=model.coef_
x=np.linspace(np.log(rsqsim.loc[condition,'SR'].to_numpy()).min(),np.log(rsqsim.loc[condition,'SR'].to_numpy()).max(),100).reshape((-1, 1))
y=model.predict(x)
print('sigma={0:.2f} x SR ^{1:.2f}'.format(np.e**a[0],b[0][0]))
axe[0].plot(np.e**x,np.e**y,'r--',linewidth=6,label='Linear Regression')
xx,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawx), statistic='mean', bins=nbin)
mean,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawy), statistic='mean', bins=nbin)
std,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawy), statistic='std', bins=nbin)
axe[0].plot(10**xx,10**mean,'-',color='darkorange',linewidth=6,label='bin-average')
sr=np.logspace(-1.5,1.6,100)
sdpm1=5.4*sr**(-0.22)
axe[0].plot(sr,sdpm1,'--',color='lightcoral',linewidth=4,label='Perrin et al (2021)')
sdpm2=30*10**(-1.17)*(sr)**(-0.226)
axe[0].plot(sr,sdpm2,'--',linewidth=4,color='gold',label='Log-scale w creeping events')
sdpm3=30*10**(-1.14)*(sr)**(-0.102)
axe[0].plot(sr,sdpm3,'--',linewidth=4,color='dodgerblue',label='Log-scale wo creeping events')
sdpm4=30*10**(-1.11)*(sr)**(-0.0769)
axe[0].plot(sr,sdpm4,'--',linewidth=4,color='green',label='Linear-scale')
axe[0].legend(loc='lower left',fontsize='x-large')
plt.tight_layout()
For reverse events¶
In [8]:
fig = plt.figure( figsize=(15,5))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='RV')
axe[0].hist(rsqsim.loc[condition,'ADip'],bins=40,density=True);
axe[1].hist(rsqsim.loc[condition,'AdDip'],bins=40,density=True);
axe[2].hist(rsqsim.loc[condition,'MdDip'],bins=40,density=True);
axe[3].hist(rsqsim.loc[condition,'ARake'],bins=40,density=True);
axe[4].hist(rsqsim.loc[condition,'AdRake'],bins=40,density=True);
axe[5].hist(rsqsim.loc[condition,'MdRake'],bins=40,density=True);
In [9]:
cmap=plt.get_cmap('GnBu',10)
sr_range=[0.01,1,10,40]
fig = plt.figure( figsize=(15,12))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='RV') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
x=rsqsim.loc[condition,'SR']
binx = np.logspace(-2,2,50)
field=['M','SDP','L','W','A','Aratio']
binyarray=[np.linspace(6.5,8.1,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[False,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('SR (mm/yr)')
plt.tight_layout()
In [10]:
cmap=plt.get_cmap('GnBu',10)
sr_range=[0.01,1,10,40]
fig = plt.figure( figsize=(15,12))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='RV') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<2))
x=rsqsim.loc[condition,'M']
binx = np.linspace(6.5,8.1,50)
field=['SR','SDP','L','W','A','Aratio']
binyarray=[np.logspace(-2,2,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[True,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
# axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('Magnitude')
plt.tight_layout()
In [11]:
cmap=plt.get_cmap('GnBu',10)
c=['darkblue','green','orange']
sr_range=[0.01,1,10,40]
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
fig = plt.figure( figsize=(15,20))
axe=fig.subplots(6, 4)
condition=(rsqsim['FM']=='RV') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
binx = np.linspace(6.5,8.3,50)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[True,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
xall=rsqsim.loc[condition,'M']
yall=rsqsim.loc[condition,field[iff]]
if logy[iff]:
yall = np.log10(yall)
meanall,bin_edges,_=stats.binned_statistic(xall, yall, statistic='mean', bins=nbin)
if logy[iff]:
meanall=10**meanall
for i in range(len(sr_range)-1):
condition2=condition & ((rsqsim['SR']>=sr_range[i]) & (rsqsim['SR']<sr_range[i+1]))
y=rsqsim.loc[condition2,field[iff]]
x=rsqsim.loc[condition2,'M']
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff,i].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
# ax1_divider = make_axes_locatable(axe[iff,i])
# cax1 = ax1_divider.append_axes("right", size="4%", pad="-10%")
# cb1 = fig.colorbar(im1, cax=cax1,orientation='vertical')
# cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff,i].set_yscale('log')
axe[iff,i].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,i].grid(which='minor', color='silver', linestyle='--')
axe[iff,i].set_axisbelow(True)
axe[iff,i].set_ylabel(field[iff],fontsize='large')
axe[iff,i].set_xlabel('Magnitude')
x = x.to_numpy()
y = y.to_numpy()
if logy[iff]:
y = np.log10(y)
if(len(y)>0):
mean,bin_edges,_=stats.binned_statistic(x, y, statistic='mean', bins=nbin)
mean_centers = (bin_edges[1:] + bin_edges[0:-1])/2
if logy[iff]:
mean=10**mean
if plotratio:
axe[iff,3].plot(mean_centers,mean/meanall,'-',color=c[i],linewidth=4,zorder=10,label='Mean for {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]) )
axe[iff,3].set_ylim([0.1,10])
axe[iff,3].set_yscale('log')
else:
if(len(y)>0):
axe[iff,3].plot(mean_centers,mean,'-',color=c[i],linewidth=4,zorder=10,label='Mean for {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]) )
if logy[iff]:
axe[iff,3].set_yscale('log')
axe[iff,3].set_ylim([biny.min(),biny.max()])
if iff==0:
axe[iff,3].legend(loc='best',fontsize='medium')
axe[iff,3].set_xlim([binx.min(),binx.max()])
axe[iff,3].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,3].grid(which='minor', color='silver', linestyle='--')
axe[iff,3].set_axisbelow(True)
plt.tight_layout()
In [12]:
cmap=plt.get_cmap('GnBu',10)
fig = plt.figure( figsize=(10,10))
axe=[]
axe.append(plt.subplot2grid((1,1), (0,0), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='RV') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<2))
x=rsqsim.loc[condition,'SR']
binx = np.logspace(-2,2,50)
field=['SDP',]
binyarray=[
np.logspace(0,1,50),
]
logy=[True,]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap,alpha=0.8,rasterized=True)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('SR (mm/yr)')
nbin=5
rawx=x.to_numpy()
rawy=y.to_numpy()
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(np.log(rsqsim.loc[condition,'SR'].to_numpy()).reshape((-1, 1)), np.log(rsqsim.loc[condition,'SDP'].to_numpy()).reshape((-1, 1)))
a=model.intercept_
b=model.coef_
x=np.linspace(np.log(rsqsim.loc[condition,'SR'].to_numpy()).min(),np.log(rsqsim.loc[condition,'SR'].to_numpy()).max(),100).reshape((-1, 1))
y=model.predict(x)
print('sigma={0:.2f} x SR ^{1:.2f}'.format(np.e**a[0],b[0][0]))
axe[0].plot(np.e**x,np.e**y,'r--',linewidth=6,label='Linear Regression')
xx,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawx), statistic='mean', bins=nbin)
mean,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawy), statistic='mean', bins=nbin)
std,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawy), statistic='std', bins=nbin)
axe[0].plot(10**xx,10**mean,'-',color='darkorange',linewidth=6,label='bin-average')
sr=np.logspace(-1.5,1.6,100)
sdpm1=5.4*sr**(-0.22)
axe[0].plot(sr,sdpm1,'--',color='lightcoral',linewidth=4,label='Perrin et al (2021)')
sdpm2=30*10**(-1.17)*(sr)**(-0.226)
axe[0].plot(sr,sdpm2,'--',linewidth=4,color='gold',label='Log-scale w creeping events')
sdpm3=30*10**(-1.14)*(sr)**(-0.102)
axe[0].plot(sr,sdpm3,'--',linewidth=4,color='dodgerblue',label='Log-scale wo creeping events')
sdpm4=30*10**(-1.11)*(sr)**(-0.0769)
axe[0].plot(sr,sdpm4,'--',linewidth=4,color='green',label='Linear-scale')
axe[0].legend(loc='lower left',fontsize='x-large')
plt.tight_layout()
Normal events¶
In [13]:
fig = plt.figure( figsize=(15,5))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='NM')
axe[0].hist(rsqsim.loc[condition,'ADip'],bins=40,density=True);
axe[1].hist(rsqsim.loc[condition,'AdDip'],bins=40,density=True);
axe[2].hist(rsqsim.loc[condition,'MdDip'],bins=40,density=True);
axe[3].hist(rsqsim.loc[condition,'ARake'],bins=40,density=True);
axe[4].hist(rsqsim.loc[condition,'AdRake'],bins=40,density=True);
axe[5].hist(rsqsim.loc[condition,'MdRake'],bins=40,density=True);
In [14]:
cmap=plt.get_cmap('Purples',10)
sr_range=[0.01,1,10,40]
fig = plt.figure( figsize=(15,12))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='NM') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
x=rsqsim.loc[condition,'SR']
binx = np.logspace(-2,2,50)
field=['M','SDP','L','W','A','Aratio']
binyarray=[np.linspace(6.5,8.1,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[False,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('SR (mm/yr)')
plt.tight_layout()
In [15]:
cmap=plt.get_cmap('Purples',10)
sr_range=[0.01,1,10,40]
fig = plt.figure( figsize=(15,12))
axe=[]
axe.append(plt.subplot2grid((2,3), (0,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,0), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,1), rowspan=1, colspan=1))
axe.append(plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='NM') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<2))
x=rsqsim.loc[condition,'M']
binx = np.linspace(6.5,8.1,50)
field=['SR','SDP','L','W','A','Aratio']
binyarray=[np.logspace(-2,2,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[True,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
# axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('Magnitude')
plt.tight_layout()
In [16]:
cmap=plt.get_cmap('Purples',10)
c=['darkblue','green','orange']
sr_range=[0.01,1,10,40]
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
fig = plt.figure( figsize=(15,20))
axe=fig.subplots(6, 4)
condition=(rsqsim['FM']=='NM') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
binx = np.linspace(6.5,8.3,50)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,50),
np.logspace(0,1,50),
np.logspace(1,3,50),
np.linspace(3,30,50),
np.logspace(2,4.1,50),
np.logspace(-1,2,50),
]
logy=[True,True,True,False,True,True]
nfield=len(field)
for iff in range(nfield):
xall=rsqsim.loc[condition,'M']
yall=rsqsim.loc[condition,field[iff]]
if logy[iff]:
yall = np.log10(yall)
meanall,bin_edges,_=stats.binned_statistic(xall, yall, statistic='mean', bins=nbin)
if logy[iff]:
meanall=10**meanall
for i in range(len(sr_range)-1):
condition2=condition & ((rsqsim['SR']>=sr_range[i]) & (rsqsim['SR']<sr_range[i+1]))
y=rsqsim.loc[condition2,field[iff]]
x=rsqsim.loc[condition2,'M']
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff,i].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap)
# ax1_divider = make_axes_locatable(axe[iff,i])
# cax1 = ax1_divider.append_axes("right", size="4%", pad="-10%")
# cb1 = fig.colorbar(im1, cax=cax1,orientation='vertical')
# cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff,i].set_yscale('log')
axe[iff,i].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,i].grid(which='minor', color='silver', linestyle='--')
axe[iff,i].set_axisbelow(True)
axe[iff,i].set_ylabel(field[iff],fontsize='large')
axe[iff,i].set_xlabel('Magnitude')
x = x.to_numpy()
y = y.to_numpy()
if logy[iff]:
y = np.log10(y)
if(len(y)>0):
mean,bin_edges,_=stats.binned_statistic(x, y, statistic='mean', bins=nbin)
mean_centers = (bin_edges[1:] + bin_edges[0:-1])/2
if logy[iff]:
mean=10**mean
if plotratio:
axe[iff,3].plot(mean_centers,mean/meanall,'-',color=c[i],linewidth=4,zorder=10,label='Mean for {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]) )
axe[iff,3].set_ylim([0.1,10])
axe[iff,3].set_yscale('log')
else:
if(len(y)>0):
axe[iff,3].plot(mean_centers,mean,'-',color=c[i],linewidth=4,zorder=10,label='Mean for {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]) )
if logy[iff]:
axe[iff,3].set_yscale('log')
axe[iff,3].set_ylim([biny.min(),biny.max()])
if iff==0:
axe[iff,3].legend(loc='best',fontsize='medium')
axe[iff,3].set_xlim([binx.min(),binx.max()])
axe[iff,3].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,3].grid(which='minor', color='silver', linestyle='--')
axe[iff,3].set_axisbelow(True)
plt.tight_layout()
In [17]:
cmap=plt.get_cmap('Purples',10)
fig = plt.figure( figsize=(10,10))
axe=[]
axe.append(plt.subplot2grid((1,1), (0,0), rowspan=1, colspan=1))
condition=(rsqsim['FM']=='NM') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<2))
x=rsqsim.loc[condition,'SR']
binx = np.logspace(-2,2,50)
field=['SDP',]
binyarray=[
np.logspace(0,1,50),
]
logy=[True,]
nfield=len(field)
for iff in range(nfield):
for i in range(len(sr_range)-1):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff].pcolormesh(xedges, yedges, np.log10(H), cmap=cmap,alpha=0.8,rasterized=True)
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("top", size="4%", pad="10%")
cb1 = fig.colorbar(im1, cax=cax1,orientation='horizontal')
cb1.set_label('$log_{10}$(density)',fontsize='large')
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_xscale('log')
axe[iff].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff].grid(which='minor', color='silver', linestyle='--')
axe[iff].set_axisbelow(True)
axe[iff].set_ylabel(field[iff],fontsize='large')
axe[iff].set_xlabel('SR (mm/yr)')
nbin=5
rawx=x.to_numpy()
rawy=y.to_numpy()
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(np.log(rsqsim.loc[condition,'SR'].to_numpy()).reshape((-1, 1)), np.log(rsqsim.loc[condition,'SDP'].to_numpy()).reshape((-1, 1)))
a=model.intercept_
b=model.coef_
x=np.linspace(np.log(rsqsim.loc[condition,'SR'].to_numpy()).min(),np.log(rsqsim.loc[condition,'SR'].to_numpy()).max(),100).reshape((-1, 1))
y=model.predict(x)
print('sigma={0:.2f} x SR ^{1:.2f}'.format(np.e**a[0],b[0][0]))
axe[0].plot(np.e**x,np.e**y,'r--',linewidth=6,label='Linear Regression')
xx,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawx), statistic='mean', bins=nbin)
mean,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawy), statistic='mean', bins=nbin)
std,_,_=stats.binned_statistic(np.log10(rawx), np.log10(rawy), statistic='std', bins=nbin)
axe[0].plot(10**xx,10**mean,'-',color='darkorange',linewidth=6,label='bin-average')
sr=np.logspace(-1.5,1.6,100)
sdpm1=5.4*sr**(-0.22)
axe[0].plot(sr,sdpm1,'--',color='lightcoral',linewidth=4,label='Perrin et al (2021)')
sdpm2=30*10**(-1.17)*(sr)**(-0.226)
axe[0].plot(sr,sdpm2,'--',linewidth=4,color='gold',label='Log-scale w creeping events')
sdpm3=30*10**(-1.14)*(sr)**(-0.102)
axe[0].plot(sr,sdpm3,'--',linewidth=4,color='dodgerblue',label='Log-scale wo creeping events')
sdpm4=30*10**(-1.11)*(sr)**(-0.0769)
axe[0].plot(sr,sdpm4,'--',linewidth=4,color='green',label='Linear-scale')
axe[0].legend(loc='lower left',fontsize='x-large')
plt.tight_layout()
In [ ]: