"""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()