probabilistic-analysis/probabilistic_assessment.py

@@ -464,12 +464,6 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                                    figsize=(16, 24),
                                    sharey='row')
 
-            # Check whether to save probabilistic diagnostics
-            for _, bp in pd.DataFrame(diagnostics).iterrows():
-                if ((str(prof['block']) == str(bp['block']))
-                        and (prof['profile'] == bp['profile'])):
-                    output_diagnostics = True
-
             # Loop through years
             pbar_year = tqdm(output_years, leave=False)
             for j, year in enumerate(pbar_year):
@@ -601,7 +595,12 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                                                  header=False,
                                                  float_format='%g')
 
-                if output_diagnostics:
+                # Check whether to save probabilistic diagnostics
+                for _, bp in pd.DataFrame(diagnostics).iterrows():
+                    if not ((str(prof['block']) == str(bp['block'])) and
+                            (prof['profile'] == bp['profile'])):
+                        continue  # Don't save
+
                 # Save probabilistic diagnostics
                 year_idx = year == years
 
@@ -658,14 +657,15 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                                               ascending=False)
 
                 # Add encounter probabilities
-                    dump_df['Encounter probability (%)'] = np.linspace(
-                        0, 100, num=n_runs + 2)[1:-1]
+                dump_df['Encounter probability (%)'] = np.linspace(0,
+                                                                   100,
+                                                                   num=n_runs +
+                                                                   2)[1:-1]
                 dump_df = dump_df.set_index('Encounter probability (%)')
 
                 csv_name = os.path.join(
                     'diagnostics',
-                        '{} {} {}.csv'.format(beach_scenario, year,
-                                              profile_type))
+                    '{} {} {}.csv'.format(beach_scenario, year, profile_type))
                 dump_df.to_csv(csv_name, float_format='%g')
 
                 for i, c in enumerate(dump_df.columns[3:]):
@@ -679,24 +679,19 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                     if j == 0:
                         ax[i, 0].yaxis.set_label_coords(-0.4, 0.5)
                         label = c.replace('(', '\n(')
-                            ax[i, 0].set_ylabel(label,
-                                                va='top',
-                                                linespacing=1.5)
+                        ax[i, 0].set_ylabel(label, va='top', linespacing=1.5)
 
-                    ax[i, j].set_xlabel('Encounter probability (%)',
-                                        labelpad=10)
+                ax[i, j].set_xlabel('Encounter probability (%)', labelpad=10)
                 ax[0, j].set_title(year)
 
                 fig.suptitle('{},  block {}, profile {}'.format(
                     beach_scenario, prof['block'], prof['profile']),
                              y=0.92)
 
-            if output_diagnostics:
             figname = os.path.join(
-                    'diagnostics',
-                    f'{beach_scenario} {profile_type} scatter.png')
+                'diagnostics', '{} {}.png'.format(beach_scenario,
+                                                  profile_type))
             plt.savefig(figname, bbox_inches='tight', dpi=300)
-                plt.close(fig)
 
 
 def main():

slr/distribution_fitting.py

@@ -1,80 +0,0 @@
-"""Fit probability distributions to IPCC sea level rise forecasts.
-
-
-Reads:
-    'IPCC AR6.xlsx'
-
-Writes:
-    'png/*.png'
-
-D. Howe
-d.howe@wrl.unsw.edu.au
-2022-05-12
-"""
-import os
-import re
-import numpy as np
-import pandas as pd
-from scipy import stats, optimize
-import matplotlib.pyplot as plt
-
-# Read data
-df = pd.read_excel('IPCC AR6.xlsx', index_col=[0, 1, 2, 3, 4])
-df = df.sort_index()
-
-# Use all 'medium' confidence scenarios for intermediate quantiles
-scenarios = ['ssp119', 'ssp126', 'ssp245', 'ssp370', 'ssp585']
-dff = df.loc[838, 'total', 'medium', scenarios].groupby('quantile').mean()
-
-# Use ssp119/ssp585 for 5th and 95th quantiles
-dff.loc[5] = df.loc[838, 'total', 'medium', 'ssp119', 5]
-dff.loc[95] = df.loc[838, 'total', 'medium', 'ssp585', 95]
-dff = dff.T
-
-dff.index.name = 'year'
-percentiles = dff.columns.values / 100
-values = dff.loc[2150].values
-
-x_min = values.min() - 0.2
-x_max = values.max() + 0.2
-x = np.linspace(x_min, x_max, num=1000)
-
-# Get statistical distributions
-distributions = [
-    getattr(stats, d) for d in dir(stats)
-    if isinstance(getattr(stats, d), stats.rv_continuous)
-]
-
-for dist in distributions:
-
-    def cdf(x, loc, scale):
-        """Calculate cumulative density function"""
-        return dist(loc=loc, scale=scale).cdf(x)
-
-    try:
-        loc, scale = optimize.curve_fit(cdf, values, percentiles)[0]
-        p = {'loc': loc, 'scale': scale}
-
-    except TypeError:
-        continue
-
-    fig, ax = plt.subplots(1, 2, figsize=(6, 2))
-
-    ax[0].plot(x, 100 * dist.cdf(x, **p))
-    ax[0].plot(values, 100 * percentiles, '.', c='#444444')
-    ax[1].plot(x, 100 * dist.pdf(x, **p))
-
-    ax[0].axhline(y=100, c='#000000', lw=0.8, zorder=-1)
-    ax[0].set_ylim(0, 101)
-    ax[1].set_ylim(bottom=0)
-
-    ax[0].set_title(dist.name, x=-0.7, y=0.5, ha='left')
-
-    for a in ax.ravel():
-        a.spines['right'].set_visible(False)
-        a.spines['top'].set_visible(False)
-
-    ax[1].set_yticks([])
-
-    plt.savefig(f'png/{dist.name}.png', bbox_inches='tight', dpi=100)
-    plt.close()