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
from sklearn.linear_model import LinearRegression
%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]:
SS=((rsqsim['FM']=='RSS') | (rsqsim['FM']=='LSS')) & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
RV=(rsqsim['FM']=='RV') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
NM=(rsqsim['FM']=='NM') & ((rsqsim['AdRake']<5) & (rsqsim['AdDip']<5))
In [3]:
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('ivory'), c('orange'),0.33,c('orange'),c('red'), 0.66,c('red'), c('black')],256)
rvb=plt.get_cmap('jet',256)
In [4]:
fm = [SS,RV,NM]
fmcmap=[rvb,
rvb,
rvb]
nx = 50
ny = 50
field=['M',
'SDP',
'L',
'W',
#'A',
'Aratio']
binx = np.logspace(-2,2,nx)
binyarray=[np.linspace(6.5,8.1,ny),
np.logspace(0,1,ny),
np.logspace(0,3,ny),
np.logspace(0.5,1.9,ny),
#np.logspace(2,4.1,50),
np.logspace(-1,2,ny)
]
logy=[False,True,True,True,True,True]
nfield=len(field)
fig = plt.figure( figsize=(15,15))
axe=fig.subplots(nfield, 3)
for ifm in range(len(fm)):
condition = fm[ifm]
x=rsqsim.loc[condition,'SR']
for iff in range(nfield):
y=rsqsim.loc[condition,field[iff]]
biny = binyarray[iff]
H, yedges, xedges = np.histogram2d(y, x, bins=[biny,binx],density=True)
im1=axe[iff,ifm].pcolormesh(xedges, yedges, np.log10(H), cmap=fmcmap[ifm])
# ax1_divider = make_axes_locatable(axe[iff,ifm])
# 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,ifm].set_yscale('log')
axe[iff,ifm].set_xscale('log')
axe[iff,ifm].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,ifm].grid(which='minor', color='silver', linestyle='--')
axe[iff,ifm].set_axisbelow(True)
axe[iff,ifm].set_ylabel(field[iff],fontsize='large')
if iff==nfield-1:
axe[iff,ifm].set_xlabel('SR (mm/yr)')
plt.tight_layout()
In [5]:
fmcmap=[rvb,
rvb,
rvb]
nm=['SS','RV','NM']
c=['red','blue','green']
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
nx = 50
ny = 50
binx = np.logspace(-2,2,nx)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,ny),
np.logspace(0,1,ny),
np.logspace(1,3,ny),
np.logspace(0.5,1.9,ny),
np.logspace(2,4.1,ny),
np.logspace(-1,2,50),
]
logy=[True,True,True,True,True,True]
nfield=len(field)
fig = plt.figure( figsize=(15,20))
axe=fig.subplots(nfield, 4)
binx = np.linspace(6.5,8.3,nx)
for iff in range(nfield):
for ifm in range(len(fm)):
condition = fm[ifm]
condition2=condition
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,ifm].pcolormesh(xedges, yedges, np.log10(H), cmap=fmcmap[ifm])
# 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,ifm].set_yscale('log')
axe[iff,ifm].grid(b=True,which='major', color='silver', linestyle='-')
axe[iff,ifm].grid(which='minor', color='silver', linestyle='--')
axe[iff,ifm].set_axisbelow(True)
axe[iff,ifm].set_ylabel(field[iff],fontsize='large')
axe[iff,ifm].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[ifm],linewidth=4,zorder=10,label=nm[ifm] )
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[ifm],linewidth=4,zorder=10,label=nm[ifm])
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()
split into different fault maturity group¶
In [6]:
fmcmap=[rvb,
rvb,
rvb]
nm=['SS','RV','NM']
c=['red','blue','green']
ls=[(0, (5, 3)),(0, (3, 2)),(0, (1, 1))]
al=[1,0.6,0.3]
zo=[10,9,8]
sr_range=[0.01,1,10,40]
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
plotratio=False
nx = 50
ny = 50
binx = np.linspace(6.5,8.3,nx)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,ny),
np.logspace(0,1,ny),
np.logspace(1,3,ny),
np.logspace(0.5,1.9,ny),
np.logspace(2,4.1,ny),
np.logspace(-1,2,50),
]
logy=[True,True,True,True,True,True]
nfield=len(field)
fig = plt.figure( figsize=(20,10))
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))
binx = np.linspace(6.5,8.3,nx)
for iff in range(nfield):
for ifm in range(len(fm)):
for i in range(len(sr_range)-1):
condition = fm[ifm]
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]
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].plot(mean_centers,mean/meanall,linestyle=ls[i],color=c[ifm],linewidth=2,zorder=10,label=nm[ifm] )
axe[iff].set_ylim([0.1,10])
axe[iff].set_yscale('log')
else:
if(len(y)>0):
axe[iff].plot(mean_centers,mean,linestyle=ls[i],color=c[ifm],linewidth=4,zorder=zo[i]-ifm,
label=nm[ifm]+' {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]),alpha=al[i])
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_ylim([biny.min(),biny.max()])
if iff==0:
axe[iff].legend(loc='best',fontsize='medium')
axe[iff].set_xlim([binx.min(),binx.max()])
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='x-large')
axe[iff].set_xlabel('Magnitude')
plt.tight_layout()
In [7]:
fmcmap=[plt.get_cmap('RdPu',256),
plt.get_cmap('GnBu',256),
plt.get_cmap('Purples',256)]
field=['SDP','SDP','SDP']
fig = plt.figure( figsize=(30,10))
axe=fig.subplots(1, len(field))
binx = np.logspace(-2,2,50)
binyarray=[
np.logspace(0,1,50),
np.logspace(0,1,50),
np.logspace(0,1,50),
]
logy=[True,True,True]
nbin=5
nfield=len(field)
for iff in range(len(fm)):
condition=fm[iff]
x=rsqsim.loc[condition,'SR']
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=fmcmap[iff])
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("bottom", 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)')
rawx=x.to_numpy()
rawy=y.to_numpy()
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)
axe[iff].set_title('sigma={0:.2f} x SR ^{1:.2f}'.format(np.e**a[0],b[0][0]))
axe[iff].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[iff].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[iff].plot(sr,sdpm1,'--',color='lightcoral',linewidth=4,label='Perrin et al (2021)')
sdpm2=30*10**(-1.17)*(sr)**(-0.226)
axe[iff].plot(sr,sdpm2,'--',linewidth=4,color='gold',label='Log-scale w creeping events')
sdpm3=30*10**(-1.14)*(sr)**(-0.102)
axe[iff].plot(sr,sdpm3,'--',linewidth=4,color='dodgerblue',label='Log-scale wo creeping events')
sdpm4=30*10**(-1.11)*(sr)**(-0.0769)
axe[iff].plot(sr,sdpm4,'--',linewidth=4,color='green',label='Linear-scale')
if iff==0:
axe[iff].legend(loc='lower left',fontsize='small')
Change defination of L¶
In [8]:
rsqsim['L']=rsqsim['FML']
rsqsim['W']=rsqsim['A']/rsqsim['L']
rsqsim['Aratio']=rsqsim['L']/rsqsim['W']
rsqsim['SF']=rsqsim['SSL']/rsqsim['L']
rsqsim['SDP']=chinnery_strdrop(rsqsim['L'].to_numpy(),rsqsim['AD'].to_numpy(),rsqsim['W'].to_numpy())
fmcmap=[plt.get_cmap('RdPu',256),
plt.get_cmap('GnBu',256),
plt.get_cmap('Purples',256)]
field=['SDP','SDP','SDP']
fig = plt.figure( figsize=(30,10))
axe=fig.subplots(1, len(field))
binx = np.logspace(-2,2,50)
binyarray=[
np.logspace(0,1,50),
np.logspace(0,1,50),
np.logspace(0,1,50),
]
logy=[True,True,True]
nbin=5
nfield=len(field)
for iff in range(len(fm)):
condition=fm[iff]
x=rsqsim.loc[condition,'SR']
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=fmcmap[iff])
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("bottom", 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)')
rawx=x.to_numpy()
rawy=y.to_numpy()
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)
axe[iff].set_title('sigma={0:.2f} x SR ^{1:.2f}'.format(np.e**a[0],b[0][0]))
axe[iff].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[iff].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[iff].plot(sr,sdpm1,'--',color='lightcoral',linewidth=4,label='Perrin et al (2021)')
sdpm2=30*10**(-1.17)*(sr)**(-0.226)
axe[iff].plot(sr,sdpm2,'--',linewidth=4,color='gold',label='Log-scale w creeping events')
sdpm3=30*10**(-1.14)*(sr)**(-0.102)
axe[iff].plot(sr,sdpm3,'--',linewidth=4,color='dodgerblue',label='Log-scale wo creeping events')
sdpm4=30*10**(-1.11)*(sr)**(-0.0769)
axe[iff].plot(sr,sdpm4,'--',linewidth=4,color='green',label='Linear-scale')
if iff==0:
axe[iff].legend(loc='lower left',fontsize='small')
In [9]:
fmcmap=[rvb,
rvb,
rvb]
nm=['SS','RV','NM']
c=['red','blue','green']
ls=[(0, (5, 3)),(0, (3, 2)),(0, (1, 1))]
al=[1,0.6,0.3]
zo=[10,9,8]
sr_range=[0.01,1,10,40]
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
plotratio=False
nx = 50
ny = 50
binx = np.linspace(6.5,8.3,nx)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,ny),
np.logspace(0,1,ny),
np.logspace(1,3,ny),
np.logspace(0.5,1.9,ny),
np.logspace(2,4.1,ny),
np.logspace(-1,2,50),
]
logy=[True,True,True,True,True,True]
nfield=len(field)
fig = plt.figure( figsize=(20,10))
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))
binx = np.linspace(6.5,8.3,nx)
for iff in range(nfield):
for ifm in range(len(fm)):
for i in range(len(sr_range)-1):
condition = fm[ifm]
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]
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].plot(mean_centers,mean/meanall,linestyle=ls[i],color=c[ifm],linewidth=2,zorder=10,label=nm[ifm] )
axe[iff].set_ylim([0.1,10])
axe[iff].set_yscale('log')
else:
if(len(y)>0):
axe[iff].plot(mean_centers,mean,linestyle=ls[i],color=c[ifm],linewidth=4,zorder=zo[i]-ifm,
label=nm[ifm]+' {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]),alpha=al[i])
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_ylim([biny.min(),biny.max()])
if iff==0:
axe[iff].legend(loc='best',fontsize='medium')
axe[iff].set_xlim([binx.min(),binx.max()])
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='x-large')
axe[iff].set_xlabel('Magnitude')
plt.tight_layout()
In [10]:
rsqsim['L']=rsqsim['MSL']
rsqsim['W']=rsqsim['A']/rsqsim['L']
rsqsim['Aratio']=rsqsim['L']/rsqsim['W']
rsqsim['SF']=rsqsim['SSL']/rsqsim['L']
rsqsim['SDP']=chinnery_strdrop(rsqsim['L'].to_numpy(),rsqsim['AD'].to_numpy(),rsqsim['W'].to_numpy())
rsqsim.dropna(inplace=True)
fmcmap=[plt.get_cmap('RdPu',256),
plt.get_cmap('GnBu',256),
plt.get_cmap('Purples',256)]
field=['SDP','SDP','SDP']
fig = plt.figure( figsize=(30,10))
axe=fig.subplots(1, len(field))
binx = np.logspace(-2,2,50)
binyarray=[
np.logspace(0,1,50),
np.logspace(0,1,50),
np.logspace(0,1,50),
]
logy=[True,True,True]
nbin=5
nfield=len(field)
for iff in range(len(fm)):
condition=fm[iff]
x=rsqsim.loc[condition,'SR']
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=fmcmap[iff])
ax1_divider = make_axes_locatable(axe[iff])
cax1 = ax1_divider.append_axes("bottom", 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)')
rawx=x.to_numpy()
rawy=y.to_numpy()
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)
axe[iff].set_title('sigma={0:.2f} x SR ^{1:.2f}'.format(np.e**a[0],b[0][0]))
axe[iff].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[iff].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[iff].plot(sr,sdpm1,'--',color='lightcoral',linewidth=4,label='Perrin et al (2021)')
sdpm2=30*10**(-1.17)*(sr)**(-0.226)
axe[iff].plot(sr,sdpm2,'--',linewidth=4,color='gold',label='Log-scale w creeping events')
sdpm3=30*10**(-1.14)*(sr)**(-0.102)
axe[iff].plot(sr,sdpm3,'--',linewidth=4,color='dodgerblue',label='Log-scale wo creeping events')
sdpm4=30*10**(-1.11)*(sr)**(-0.0769)
axe[iff].plot(sr,sdpm4,'--',linewidth=4,color='green',label='Linear-scale')
if iff==0:
axe[iff].legend(loc='lower left',fontsize='small')
In [11]:
fmcmap=[rvb,
rvb,
rvb]
nm=['SS','RV','NM']
c=['red','blue','green']
ls=[(0, (5, 3)),(0, (3, 2)),(0, (1, 1))]
al=[1,0.6,0.3]
zo=[10,9,8]
sr_range=[0.01,1,10,40]
nbin=np.arange(6.5,8.2,0.2)
plotratio=False
plotratio=False
nx = 50
ny = 50
binx = np.linspace(6.5,8.3,nx)
field=['AD','SDP','L','W','A','Aratio']
binyarray=[np.logspace(0,1,ny),
np.logspace(0,1,ny),
np.logspace(1,3,ny),
np.logspace(0.5,1.9,ny),
np.logspace(2,4.1,ny),
np.logspace(-1,2,50),
]
logy=[True,True,True,True,True,True]
nfield=len(field)
fig = plt.figure( figsize=(20,10))
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))
binx = np.linspace(6.5,8.3,nx)
for iff in range(nfield):
for ifm in range(len(fm)):
for i in range(len(sr_range)-1):
condition = fm[ifm]
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]
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].plot(mean_centers,mean/meanall,linestyle=ls[i],color=c[ifm],linewidth=2,zorder=10,label=nm[ifm] )
axe[iff].set_ylim([0.1,10])
axe[iff].set_yscale('log')
else:
if(len(y)>0):
axe[iff].plot(mean_centers,mean,linestyle=ls[i],color=c[ifm],linewidth=4,zorder=zo[i]-ifm,
label=nm[ifm]+' {0}<SR<{1}'.format(sr_range[i],sr_range[i+1]),alpha=al[i])
if logy[iff]:
axe[iff].set_yscale('log')
axe[iff].set_ylim([biny.min(),biny.max()])
if iff==0:
axe[iff].legend(loc='best',fontsize='medium')
axe[iff].set_xlim([binx.min(),binx.max()])
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='x-large')
axe[iff].set_xlabel('Magnitude')
plt.tight_layout()
In [ ]: