From f66e83be64f19787e1c3466908e1a9271297e95a Mon Sep 17 00:00:00 2001 From: Dan Howe Date: Sun, 15 May 2022 17:13:30 +0200 Subject: [PATCH] Plot diagnostics time series --- .../probabilistic_assessment.py | 121 +++++++++++++++++- 1 file changed, 116 insertions(+), 5 deletions(-) diff --git a/probabilistic-analysis/probabilistic_assessment.py b/probabilistic-analysis/probabilistic_assessment.py index 98a97d1..903f164 100644 --- a/probabilistic-analysis/probabilistic_assessment.py +++ b/probabilistic-analysis/probabilistic_assessment.py @@ -475,17 +475,18 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year, xp=profile_volume[valid_idx], fp=profile_chainage[valid_idx]) - fig, ax = plt.subplots(9, - len(output_years), - figsize=(16, 24), - sharey='row') - # Check whether to save probabilistic diagnostics + output_diagnostics = False for _, bp in pd.DataFrame(diagnostics).iterrows(): if ((str(prof['block']) == str(bp['block'])) and (prof['profile'] == bp['profile'])): output_diagnostics = True + fig, ax = plt.subplots(9, + len(output_years), + figsize=(16, 24), + sharey='row') + # Loop through years pbar_year = tqdm(output_years, leave=False) for j, year in enumerate(pbar_year): @@ -682,6 +683,7 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year, 'diagnostics', '{} {} {}.csv'.format(beach_scenario, year, profile_type)) + dump_df = dump_df[::100] # Only output every 100th row dump_df.to_csv(csv_name, float_format='%g') for i, c in enumerate(dump_df.columns[3:]): @@ -714,6 +716,115 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year, plt.savefig(figname, bbox_inches='tight', dpi=300) plt.close(fig) + # Plot time series figure + fig, ax = plt.subplots(4, + 2, + figsize=(12, 16), + sharey='row', + gridspec_kw={ + 'wspace': 0.05, + 'width_ratios': [3, 1] + }) + + ax[0, 0].plot(years, slr[:, :100], c='#aaaaaa', lw=0.2) + ax[0, 0].plot(years, + slr[:, 1], + c='C0', + label='Sample simulation') + ax[0, 1].hist(slr[-1, :], + bins=100, + fc='#cccccc', + ec='#aaaaaa', + orientation='horizontal') + + ax[1, 0].plot(years, (slr * bf)[:, :100], c='#aaaaaa', lw=0.2) + ax[1, 0].plot(years, (slr * bf)[:, 1], c='C0') + ax[1, 1].hist((slr * bf)[-1, :], + bins=100, + fc='#cccccc', + ec='#aaaaaa', + orientation='horizontal') + + ax[2, 0].plot(years, ur[:, :100], c='#aaaaaa', lw=0.2) + ax[2, 0].plot(years, ur[:, 1], c='C0') + ax[2, 1].hist(ur[-1, :], + bins=100, + fc='#cccccc', + ec='#aaaaaa', + orientation='horizontal') + + ax[3, 0].plot(years, r[:, :100], c='#aaaaaa', lw=0.2) + ax[3, 0].plot(years, r[:, 1], c='C0', zorder=3) + ax[3, 1].hist(r[-1, :], + bins=100, + fc='#cccccc', + ec='#aaaaaa', + orientation='horizontal') + + baseline = r[:, 1] + sdd = storm_demand_dist[:, 1] + for i in range(len(slr)): + pe = [ + matplotlib.patheffects.Stroke(linewidth=5, + foreground='#ffffff', + capstyle='butt'), + matplotlib.patheffects.Normal() + ] + + ax[3, 0].plot([years[i], years[i]], + [baseline[i], baseline[i] + sdd[i]], + c='C0', + path_effects=pe) + + # Maximum recession encountered + r_max = (baseline + sdd).max() + ax[3, 0].axhline(y=r_max, c='C3', linestyle=':') + + i = len(years) - 1 + for a in ax[:, 0]: + a.set_xlim(right=years[-1]) + + # Add line at zero + for a in ax[:-1, 0]: + a.spines['bottom'].set_position(('data', 0)) + a.set_xticklabels([]) + + for a in ax[:, 1]: + a.xaxis.set_visible(False) + a.spines['bottom'].set_visible(False) + + ax[3, 0].annotate( + 'Most eroded beach state encountered in planning period', + (years[0], r_max), + xytext=(0, 20), + textcoords='offset pixels') + + ax[0, 0].legend() + + ax[0, 0].set_title((f'Probabilistic trajectories\n' + f'(first 100 out of {n_runs:,} runs)')) + ax[0, 1].set_title( + f'Probability\ndistribution\nin year {years[i]}') + ax[0, 0].set_ylabel('Sea level (m)', labelpad=10) + ax[1, 0].set_ylabel('Bruun recession (m)', labelpad=10) + ax[2, 0].set_ylabel('Underlying recession (m)', labelpad=10) + ax[3, 0].set_ylabel('Shoreline displacement (m)', labelpad=10) + + ax[3, 0].set_title(('Bruun recession' + '+ underlying recession' + '+ storm demand'), + y=0.9) + + for a in ax.ravel(): + a.spines['top'].set_visible(False) + a.spines['right'].set_visible(False) + + figname = os.path.join( + 'diagnostics', + f'{beach_scenario} {profile_type} timeseries.png') + plt.savefig(figname, bbox_inches='tight', dpi=300) + plt.close(fig) + def main():