import sys
sys.path.append('../eq_stats/')
import numpy as np
import pandas as pd
import eq_stats as eqs
Define some parameters for calculations, and import fault data
M_vec = np.linspace(5, 7.64, num=1000)
FM_vec = eqs.F(M_vec)
n_samp = 1e5
f = pd.read_csv('../data/lanf_stats.csv', index_col=0)
f
L_km | z_km | slip_rate_mm_a | sr_err_mm_a | dip_deg | dip_err_deg | |
---|---|---|---|---|---|---|
fault | ||||||
kongur_shan | 147 | 10 | 6.5 | 2.0 | 30 | 10.0 |
leo_pargil | 56 | 10 | 1.0 | 1.0 | 25 | 10.0 |
gurla_mandhata | 44 | 10 | 4.5 | 2.0 | 33 | 10.0 |
s_lunggar | 20 | 10 | 2.5 | 0.5 | 22 | 5.0 |
n_lunggar | 32 | 10 | 2.0 | 1.0 | 25 | 5.0 |
pqx_qingdu | 15 | 10 | 1.0 | 0.5 | 30 | 5.0 |
pqx_n | 11 | 10 | 1.0 | 0.5 | 25 | 7.5 |
nqtl | 115 | 10 | 2.5 | 1.5 | 30 | 8.0 |
alto_tiberina | 60 | 4 | 2.4 | 0.3 | 25 | 5.0 |
time_window = np.arange(100) + 1
mag_vals = np.array([5.0, 5.5, 6.0, 6.5, 7.0, 7.5])
make fault dataframe
slr = f.loc['s_lunggar']
df_cols = ['M', 'Mo', 'Ddot', 'dip', 'recur_int', 'cum_yrs']
slr_mc = pd.DataFrame(index=np.arange(n_samp), columns=df_cols)
Populate dataframe with sample of eqs, dips, slip rates
slr_mc.M = eqs.sample_from_pdf(M_vec, FM_vec, n_samp)
slr_mc.Mo = eqs.calc_Mo_from_M(slr_mc.M)
slr_mc.dip, prob_frac = eqs.dip_rand_samp(slr['dip_deg'], slr['dip_err_deg'],
n_samp)
slr_mc.Ddot = eqs.Ddot_rand_samp(slr['slip_rate_mm_a'], slr['sr_err_mm_a'],
n_samp)
slr_mc.recur_int = eqs.calc_recurrence_interval(Mo=slr_mc.Mo, dip=slr_mc.dip,
slip_rate=slr_mc.Ddot,
L=slr['L_km'], z=slr['z_km'])
slr_mc.cum_yrs = eqs.calc_cumulative_yrs(slr_mc.recur_int)
Make EQ time series and calculate probability
slr_eq_time_series = eqs.make_eq_time_series(slr_mc.M, slr_mc.cum_yrs)
eqs.get_prob_above_val_in_window(slr_eq_time_series, 6, 40)
0.023580344543560616
make dataframe of time windows, magnitudes, and probabilities
slr_probs = pd.DataFrame(index=time_window, columns=mag_vals)
for t in time_window:
for M in mag_vals:
#print t, M
slr_probs.loc[t][M] = eqs.get_prob_above_val_in_window(
slr_eq_time_series, M, t) * prob_frac
slr_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.005726721 | 0.001839079 | 0.0005895086 | 0.0001850303 | 4.9078e-05 | 5.669453e-06 |
2 | 0.01145338 | 0.003678158 | 0.001179017 | 0.0003700607 | 9.815599e-05 | 1.133891e-05 |
3 | 0.01718005 | 0.005517237 | 0.001768526 | 0.000555091 | 0.000147234 | 1.700836e-05 |
4 | 0.02280644 | 0.007356316 | 0.002358034 | 0.0007401214 | 0.000196312 | 2.267781e-05 |
5 | 0.02798162 | 0.009195395 | 0.002947543 | 0.0009251517 | 0.00024539 | 2.834727e-05 |
!mkdir ../results
slr_probs.to_csv('../results/slr_table.csv')
def get_probs(fault_name):
#pr = prefix
# get fault info
pr = f.loc[fault_name]
df_cols = ['M', 'Mo', 'Ddot', 'dip', 'recur_int', 'cum_yrs']
pr_mc = pd.DataFrame(index=np.arange(n_samp), columns=df_cols)
# Populate dataframe with sample of eqs, dips, slip rates
pr_mc.M = eqs.sample_from_pdf(M_vec, FM_vec, n_samp)
pr_mc.Mo = eqs.calc_Mo_from_M(pr_mc.M)
pr_mc.dip, prob_frac = eqs.dip_rand_samp(pr['dip_deg'], pr['dip_err_deg'],
n_samp)
pr_mc.Ddot = eqs.Ddot_rand_samp(pr['slip_rate_mm_a'], pr['sr_err_mm_a'],
n_samp)
pr_mc.recur_int = eqs.calc_recurrence_interval(Mo=pr_mc.Mo, dip=pr_mc.dip,
slip_rate=pr_mc.Ddot,
L=pr['L_km'], z=pr['z_km'])
pr_mc.cum_yrs = eqs.calc_cumulative_yrs(pr_mc.recur_int)
# Make EQ time series and calculate probability
pr_eq_time_series = eqs.make_eq_time_series(pr_mc.M, pr_mc.cum_yrs)
# make dataframe of time windows, magnitudes, and probabilities
pr_probs = pd.DataFrame(index=time_window, columns=mag_vals)
for t in time_window:
for M in mag_vals:
pr_probs.loc[t][M] = eqs.get_prob_above_val_in_window(
pr_eq_time_series, M, t) * prob_frac
return pr_probs
nlr_probs = get_probs('n_lunggar')
nlr_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.005830359 | 0.001880466 | 0.0006014599 | 0.0001904195 | 5.078243e-05 | 6.355091e-06 |
2 | 0.01166066 | 0.003760873 | 0.001202861 | 0.0003807808 | 0.0001015649 | 1.271018e-05 |
3 | 0.01744269 | 0.005641281 | 0.001804263 | 0.000571142 | 0.0001523473 | 1.906527e-05 |
4 | 0.0228444 | 0.007521688 | 0.002405664 | 0.0007615032 | 0.0002031297 | 2.542037e-05 |
5 | 0.02775362 | 0.009402096 | 0.003007066 | 0.0009518644 | 0.0002539121 | 3.177546e-05 |
nlr_probs.to_csv('../results/nlr_table.csv')
ks_probs = get_probs('kongur_shan')
ks_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.02936775 | 0.009361432 | 0.002965765 | 0.0009282716 | 0.0002471871 | 3.023842e-05 |
2 | 0.0436878 | 0.0184368 | 0.005931531 | 0.001856543 | 0.0004943742 | 6.047685e-05 |
3 | 0.05386519 | 0.02665253 | 0.008897296 | 0.002784815 | 0.0007415613 | 9.071527e-05 |
4 | 0.06201085 | 0.03404078 | 0.01186306 | 0.003713086 | 0.0009887484 | 0.0001209537 |
5 | 0.06894744 | 0.04072156 | 0.01482883 | 0.004641358 | 0.001235936 | 0.0001511921 |
ks_probs.to_csv('../results/ks_table.csv')
lp_probs = get_probs('leo_pargil')
lp_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.0008126479 | 0.0002593972 | 8.3077e-05 | 2.629729e-05 | 7.053784e-06 | 7.720155e-07 |
2 | 0.001624361 | 0.0005187944 | 0.000166154 | 5.259457e-05 | 1.410757e-05 | 1.544031e-06 |
3 | 0.002404544 | 0.0007781916 | 0.000249231 | 7.889186e-05 | 2.116135e-05 | 2.316047e-06 |
4 | 0.003127191 | 0.001037589 | 0.000332308 | 0.0001051891 | 2.821514e-05 | 3.088062e-06 |
5 | 0.003795813 | 0.001296986 | 0.000415385 | 0.0001314864 | 3.526892e-05 | 3.860078e-06 |
lp_probs.to_csv('../results/lp_table.csv')
gm_probs = get_probs('gurla_mandhata')
gm_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.006218659 | 0.002000729 | 0.0006392782 | 0.0001975046 | 5.360484e-05 | 6.716152e-06 |
2 | 0.0115432 | 0.004001458 | 0.001278556 | 0.0003950092 | 0.0001072097 | 1.34323e-05 |
3 | 0.01549939 | 0.006002188 | 0.001917834 | 0.0005925138 | 0.0001608145 | 2.014846e-05 |
4 | 0.01867637 | 0.008002917 | 0.002557113 | 0.0007900184 | 0.0002144194 | 2.686461e-05 |
5 | 0.02136986 | 0.0100034 | 0.003196391 | 0.0009875231 | 0.0002680242 | 3.358076e-05 |
gm_probs.to_csv('../results/gm_table.csv')
pqxq_probs = get_probs('pqx_qingdu')
pqxq_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.00064289 | 0.000207872 | 6.751631e-05 | 2.054676e-05 | 5.567427e-06 | 6.236033e-07 |
2 | 0.001285774 | 0.0004157441 | 0.0001350326 | 4.109353e-05 | 1.113485e-05 | 1.247207e-06 |
3 | 0.001928657 | 0.0006236161 | 0.0002025489 | 6.164029e-05 | 1.670228e-05 | 1.87081e-06 |
4 | 0.002571541 | 0.0008314882 | 0.0002700652 | 8.218706e-05 | 2.226971e-05 | 2.494413e-06 |
5 | 0.003214424 | 0.00103936 | 0.0003375815 | 0.0001027338 | 2.783714e-05 | 3.118016e-06 |
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.0006412814 | 0.0002053575 | 6.616742e-05 | 2.063002e-05 | 5.694579e-06 | 6.028045e-07 |
2 | 0.001282556 | 0.0004107151 | 0.0001323348 | 4.126005e-05 | 1.138916e-05 | 1.205609e-06 |
3 | 0.001923831 | 0.0006160726 | 0.0001985022 | 6.189007e-05 | 1.708374e-05 | 1.808414e-06 |
4 | 0.002565106 | 0.0008214302 | 0.0002646697 | 8.252009e-05 | 2.277832e-05 | 2.411218e-06 |
5 | 0.003206381 | 0.001026788 | 0.0003308371 | 0.0001031501 | 2.847289e-05 | 3.014023e-06 |
pqxq_probs.to_csv('../results/pqxq_table.csv')
pqxn_probs = get_probs('pqx_n')
pqxn_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.0008962962 | 0.0002869313 | 9.040043e-05 | 2.783e-05 | 7.860517e-06 | 9.590369e-07 |
2 | 0.001792583 | 0.0005738626 | 0.0001808009 | 5.565999e-05 | 1.572103e-05 | 1.918074e-06 |
3 | 0.002688871 | 0.0008607939 | 0.0002712013 | 8.348999e-05 | 2.358155e-05 | 2.877111e-06 |
4 | 0.003585158 | 0.001147725 | 0.0003616017 | 0.00011132 | 3.144207e-05 | 3.836148e-06 |
5 | 0.004481445 | 0.001434656 | 0.0004520022 | 0.00013915 | 3.930259e-05 | 4.795184e-06 |
pqxn_probs.to_csv('../results/pqxn_table.csv')
nqtl_probs = get_probs('nqtl')
nqtl_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.01217447 | 0.003928182 | 0.001240889 | 0.0003806581 | 0.0001058671 | 1.378224e-05 |
2 | 0.02072983 | 0.007856242 | 0.002481657 | 0.0007611943 | 0.0002117342 | 2.756448e-05 |
3 | 0.02695721 | 0.0117843 | 0.003722424 | 0.00114173 | 0.0003176013 | 4.134672e-05 |
4 | 0.03197577 | 0.01568846 | 0.004963192 | 0.001522267 | 0.0004234684 | 5.512896e-05 |
5 | 0.03623558 | 0.01951821 | 0.006203959 | 0.001902803 | 0.0005293356 | 6.89112e-05 |
nqtl_probs.to_csv('../results/nqtl_table.csv')
at_probs = get_probs('alto_tiberina')
at_probs.head()
5.0 | 5.5 | 6.0 | 6.5 | 7.0 | 7.5 | |
---|---|---|---|---|---|---|
1 | 0.005915221 | 0.001912036 | 0.0006129352 | 0.0001919489 | 4.93921e-05 | 6.033526e-06 |
2 | 0.01183038 | 0.003824072 | 0.00122587 | 0.0003838979 | 9.87842e-05 | 1.206705e-05 |
3 | 0.01774555 | 0.005736108 | 0.001838806 | 0.0005758468 | 0.0001481763 | 1.810058e-05 |
4 | 0.02358209 | 0.007648145 | 0.002451741 | 0.0007677957 | 0.0001975684 | 2.41341e-05 |
5 | 0.02889586 | 0.009560181 | 0.003064676 | 0.0009597447 | 0.0002469605 | 3.016763e-05 |
at_probs.to_csv('../results/at_table.csv')