Commit d986d221 authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

Modified plot k func

parent 92362d88
...@@ -163,7 +163,7 @@ def filter_10yr_01_2010(cat): ...@@ -163,7 +163,7 @@ def filter_10yr_01_2010(cat):
start_date = time_utils.strptime_to_utc_epoch('2010-01-01 00:00:00.00') start_date = time_utils.strptime_to_utc_epoch('2010-01-01 00:00:00.00')
end_date = time_utils.strptime_to_utc_epoch('2020-01-01 00:00:00.00') end_date = time_utils.strptime_to_utc_epoch('2020-01-01 00:00:00.00')
min_mag = 4.95 min_mag = 4.95
max_depth = 500 max_depth = 30
cat.filter([f'origin_time >= {start_date}', cat.filter([f'origin_time >= {start_date}',
f'origin_time < {end_date}', f'origin_time < {end_date}',
......
...@@ -43,7 +43,7 @@ def K_r(Events, rates, Dists, r): ...@@ -43,7 +43,7 @@ def K_r(Events, rates, Dists, r):
#################### CATALOG #################### CATALOG
plt.close('all') plt.close('all')
cat = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino()) cat = get_catalogs.filter_10yr_01_2010(get_catalogs.emrcmt())
region = cat.region region = cat.region
counts= cat.spatial_counts() counts= cat.spatial_counts()
...@@ -59,11 +59,10 @@ Events = np.array([[i, j] for i,j in zip(cat.get_longitudes(), cat.get_latitudes ...@@ -59,11 +59,10 @@ Events = np.array([[i, j] for i,j in zip(cat.get_longitudes(), cat.get_latitudes
## Forecast ## Forecast
models = get_models.ten_years() models = get_models.ten_years()
R_s = np.arange(1, 400, 10) R_s = np.arange(0, 400, 20)
for model_i in range(19): for model_i in range(6):
plt.figure() plt.figure(figsize=(6,4))
# model_i = 11 plt.title('$K-$function\n %s' % models[model_i].name)
plt.title(models[model_i].name)
forecast_data = models[model_i].spatial_counts() forecast_data = models[model_i].spatial_counts()
rates_eqk = forecast_data[cat_ind] rates_eqk = forecast_data[cat_ind]
Dists = dists(Events) Dists = dists(Events)
...@@ -81,9 +80,10 @@ for model_i in range(19): ...@@ -81,9 +80,10 @@ for model_i in range(19):
scale = n_obs / n_fore scale = n_obs / n_fore
expected_forecast_count = int(n_obs) expected_forecast_count = int(n_obs)
num_events_to_simulate = expected_forecast_count num_events_to_simulate = expected_forecast_count
print(num_events_to_simulate)
sampling_weights = numpy.cumsum(forecast_data.ravel()) / numpy.sum(forecast_data) sampling_weights = numpy.cumsum(forecast_data.ravel()) / numpy.sum(forecast_data)
n_sim = 100 n_sim = 50
K_sim = [] K_sim = []
for i in range(n_sim): for i in range(n_sim):
...@@ -96,7 +96,6 @@ for model_i in range(19): ...@@ -96,7 +96,6 @@ for model_i in range(19):
K = [] K = []
for i in R_s: for i in R_s:
K.append(K_r(events_sim, rates_sim, dists_sim, i)) K.append(K_r(events_sim, rates_sim, dists_sim, i))
# plt.plot(R_s, K, alpha=0.1, color='red')
K_sim.append(K) K_sim.append(K)
K_sim = np.array(K_sim) K_sim = np.array(K_sim)
...@@ -115,15 +114,17 @@ for model_i in range(19): ...@@ -115,15 +114,17 @@ for model_i in range(19):
L = np.divide(np.sqrt(K_sim), K_std) L = np.divide(np.sqrt(K_sim), K_std)
L_avg = np.mean(L, axis=0) L_avg = np.mean(L, axis=0)
L_up = np.quantile(L, 0.95, axis=0) L_up = np.quantile(L, 0.9, axis=0)
L_down = np.quantile(L, 0.05, axis=0) L_down = np.quantile(L, 0.1, axis=0)
L_cat = np.divide(np.sqrt(K_cat), K_std) L_cat = np.divide(np.sqrt(K_cat), K_std)
plt.plot(R_s, L_cat, color='red') plt.plot(R_s, L_cat, color='red')
plt.plot(R_s, L_avg) plt.plot(R_s, L_avg)
plt.fill_between(R_s, L_down, L_up, plt.fill_between(R_s, L_down, L_up,
color='gray', alpha=0.2) color='gray', alpha=0.2)
plt.xlabel('$r$ [km]')
plt.ylabel('$L(r)=\frac{K(r)}{\pi}$')
plt.grid() plt.grid()
plt.show() plt.show()
# cmap = plt.get_cmap('jet') # cmap = plt.get_cmap('jet')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment