diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index b742904..d498b66 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -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)