official-apps/SeismicHazardForecasting
10
0
Fork 1
Code Issues 1 Pull Requests Activity

Update src/seismic_hazard_forecasting.py

Browse Source
This commit is contained in:
ftong 2025-06-30 18:13:04 +02:00
parent 1d6002e818
commit 3d12388da9
1 changed files with 109 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
src/seismic_hazard_forecasting.py
View File

@ -1,5 +1,113 @@
# -*- coding: utf-8 -*-
def calc_bins2(dates, unit, bin_dur, dates_calc, rate_forecast, rate_unc_high, rate_unc_low, multiplicator, filename = "activity_rate.png", **kwargs):
start_date = dates.min()
end_date = dates.max()
end_date_dt=datenum_to_datetime(end_date)
# Unit and duration check, default to 15 DAYS or 12 HOURS
# valid_units = ['hours', 'days']
# if unit not in valid_units:
# unit = 'days'
# bin_dur = 15
# Adjust to '12 hours' if data duration is less than 15 days
# if (end_date - start_date) < 15 and unit == 'days':
# unit = 'hours'
# bin_dur = 12
# -------- bin, bin edges, bin midpoint calc -------------
# Determine number of bins
bin_edges=[end_date]
while bin_edges[-1] > start_date:
bin_edges.append(bin_edges[-1] - (bin_dur/multiplicator))
bin_edges = bin_edges[::-1]
bin_edges_dt = [datenum_to_datetime(d) for d in bin_edges]
# Discretize the time data
bin_counts, _ = np.histogram(dates, bins=bin_edges)
# Calculate the rate
act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator) * multiplicator for i, count in enumerate(bin_counts)]
# ---------Apply BEAST (Zhao et al., 2019) 10 times--------
ppr_all = []
for _ in range(10):
out, ppr = apply_beast(act_rate, **kwargs) # Apply BEAST
ppr_iter = np.concatenate([p for p in ppr if p is not None], axis=0) # Concatenate results
ppr_all.append(ppr_iter[:, 1])
# Concatenate all probabilities from all BEAST runs
pprs = np.concatenate(ppr_all, axis=0)
pprs = pprs[~np.isnan(pprs)]
# -------Changepoint decision and rate calculations--------
rt, idx, u_e = av_rates(act_rate, pprs)
# --------------------- Line Plot -------------------------
fig, ax = plt.subplots(figsize=(14, 5))
# Plot activiry rate for all catalog
ax.plot(bin_edges_dt[1:len(bin_edges)], act_rate, '-o', linewidth=2, markersize=6)
# Plot forecasted value with uncertainties with red line connected to last value
next_date = end_date + (bin_dur/multiplicator)
ax.plot(datenum_to_datetime(next_date), rate_forecast, 'ro', label='Forecasted Rate', markersize=6)
ax.plot([bin_edges_dt[-1], datenum_to_datetime(next_date)], [act_rate[-1], rate_forecast], 'r-')
ax.vlines(datenum_to_datetime(next_date), rate_unc_low, rate_unc_high, colors='r', linewidth=2, label='Bootstrap uncertainty') # Vertical line for the uncertainty range
# Format the x-axis with dates
ax.xaxis.set_major_locator(AutoDateLocator())
ax.xaxis.set_major_formatter(DateFormatter('%d-%b-%Y'))
plt.xticks(rotation=45)
plt.title(f'Activity rate (Time Unit: {unit}, Time Window Duration: {bin_dur} {unit})')
plt.xlabel('Time (Window Center Date)')
plt.ylabel('Activity rate per selected time window')
plt.grid(True)
# Add Vertical Line and Shaded Area
sel_data = bin_edges[-3]
sel_data_dt=datenum_to_datetime(sel_data)
ax.axvline(sel_data_dt, color='r', linestyle='--', linewidth=2, label='Time window used')
def convert_dates_to_num(dates):
return [mdates.date2num(date) if isinstance(date, datetime) else date for date in dates]
# Shade the area to the right of the vertical line
ylim = ax.get_ylim()
# Add the shaded area using the Polygon function
ax.fill_between([sel_data_dt, end_date_dt], 0.0, ylim[1]+0.1, color='r', alpha=0.5)
# Add rates and changepoints from BEAST - changepoint belongs to NEXT period
if 'rt' in locals():
for i in range(len(idx)):
# Plot first period
if i == 0 and len(idx) > 1:
ax.plot(bin_edges_dt[1:len(bin_edges)][:int(u_e[idx[i]])], [rt[i]] * int(u_e[idx[i]]), linewidth=2)
# Plot Last and second to last periods when more than 1 changepoints
elif i == len(idx) - 1 and len(idx) > 1:
ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i - 1]]):int(u_e[idx[i]])],
[rt[i]] * (int(u_e[idx[i]]) - int(u_e[idx[i - 1]])), linewidth=2)
ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i]]):],
[rt[i + 1]] * (len(act_rate) - int(u_e[idx[i]])), linewidth=2)
# Plot first and last period if only one changepoint
elif i == 0 and len(idx) == 1:
ax.plot(bin_edges_dt[1:len(bin_edges)][:int(u_e[idx[i]])], [rt[i]] * int(u_e[idx[i]]), linewidth=2)
ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i]]):],
[rt[i + 1]] * (len(act_rate) - int(u_e[idx[i]])), linewidth=2)
# Plot intermediate value if it's only one
elif len(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[1 - 1]]):int(u_e[idx[1]])])==1:
nnn = multiplicator / bin_dur
ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i - 1]]):int(u_e[idx[i]]) + 1],
[rt[i]] * (int(u_e[idx[i]]) + 1 - int(u_e[idx[i - 1]])), linewidth=2)
# Plot intermediate periods
else:
ax.plot(bin_edges_dt[1:len(bin_edges)][int(u_e[idx[i - 1]]):int(u_e[idx[i]])],
[rt[i]] * (int(u_e[idx[i]]) - int(u_e[idx[i - 1]])), linewidth=2)
plt.legend(loc='best')
plt.savefig(filename, dpi=600)
plt.show()
return act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e
def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max,
m_kde_method, xy_select, grid_dim, xy_win_method, rate_select, time_win_duration,
@ -361,9 +469,7 @@ verbose: {verbose}")
rate_forecast = multiplicator / np.median(forecasts) # [per time unit]
# Plot of forecasted activity rate with previous binned activity rate
np.save('KDE_magnitude_PDF.csv', np.array(datenum_data))
sys.exit(0)
act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins(np.array(datenum_data), time_unit,
act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins2(np.array(datenum_data), time_unit,
time_win_duration, dates_calc,
rate_forecast, rate_unc_high, rate_unc_low,
multiplicator, quiet=True)