danh
/
roches-probabilistic-hazard-lines
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Compare commits

base: danh:18d7f448a483d474d90307a894503e00b542fabc
Branches Tags
danh:migrationtestbranch
danh:master
..
compare: danh:bf916f5f25090f8b1d44c141e62b7ab57e56e5f7
Branches Tags
danh:migrationtestbranch
danh:master

No commits in common. '18d7f448a483d474d90307a894503e00b542fabc' and 'bf916f5f25090f8b1d44c141e62b7ab57e56e5f7' have entirely different histories.
18d7f448a4 ... bf916f5f25

4 changed files with 96 additions and 181 deletions
Show Stats

0
slr/IPCC AR6.xlsx → inputs/IPCC AR6.xlsx
Unescape Escape View File

0
slr/distributions.py → inputs/distributions.py
Unescape Escape View File

197
probabilistic-analysis/probabilistic_assessment.py
Unescape Escape View File

@ -464,12 +464,6 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
figsize=(16, 24), figsize=(16, 24),
sharey='row') 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 # Loop through years
pbar_year = tqdm(output_years, leave=False) pbar_year = tqdm(output_years, leave=False)
for j, year in enumerate(pbar_year): for j, year in enumerate(pbar_year):
@ -601,102 +595,103 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
header=False, header=False,
float_format='%g') float_format='%g')
if output_diagnostics: # Check whether to save probabilistic diagnostics
# Save probabilistic diagnostics for _, bp in pd.DataFrame(diagnostics).iterrows():
year_idx = year == years if not ((str(prof['block']) == str(bp['block'])) and
(prof['profile'] == bp['profile'])):
# Find index where most extreme event occurred continue # Don't save
event_year_idx = chainage_with_recession.argmin(axis=0)
# Save probabilistic diagnostics
# define dummy index year_idx = year == years
ix = np.arange(n_runs)
dump_data = { # Find index where most extreme event occurred
'Sea level rise (m)': event_year_idx = chainage_with_recession.argmin(axis=0)
slr[event_year_idx, ix].ravel(),
'Bruun factor (-)': # define dummy index
bf[event_year_idx, ix].ravel(), ix = np.arange(n_runs)
'Bruun factor x SLR (m)': dump_data = {
slr[event_year_idx, ix].ravel() * 'Sea level rise (m)':
bf[event_year_idx, ix].ravel(), slr[event_year_idx, ix].ravel(),
'Underlying trend rate (m/yr)': 'Bruun factor (-)':
ur_rate[year_idx, :].ravel(), bf[event_year_idx, ix].ravel(),
'Underlying trend (m)': 'Bruun factor x SLR (m)':
ur[event_year_idx, ix].ravel(), slr[event_year_idx, ix].ravel() *
'Underlying + SLR (m)': bf[event_year_idx, ix].ravel(),
r[event_year_idx, ix].ravel(), 'Underlying trend rate (m/yr)':
'Total movement (m)': ur_rate[year_idx, :].ravel(),
(storm_demand_dist + r)[event_year_idx, ix].ravel(), 'Underlying trend (m)':
'Storm demand distance (m)': ur[event_year_idx, ix].ravel(),
storm_demand_dist[event_year_idx, ix].ravel(), 'Underlying + SLR (m)':
'Storm demand volume (m3/m)': r[event_year_idx, ix].ravel(),
storm_demand_volume[event_year_idx, ix].ravel(), 'Total movement (m)':
} (storm_demand_dist + r)[event_year_idx, ix].ravel(),
'Storm demand distance (m)':
dump_df = pd.DataFrame(dump_data) storm_demand_dist[event_year_idx, ix].ravel(),
dump_df['Run ID'] = np.arange(len(event_year_idx)) + 1 'Storm demand volume (m3/m)':
dump_df['Event year'] = years[event_year_idx] storm_demand_volume[event_year_idx, ix].ravel(),
dump_df['Years elapsed'] = event_year_idx + 1 }
# Reorder columns dump_df = pd.DataFrame(dump_data)
dump_df = dump_df[[ dump_df['Run ID'] = np.arange(len(event_year_idx)) + 1
'Run ID', dump_df['Event year'] = years[event_year_idx]
'Event year', dump_df['Years elapsed'] = event_year_idx + 1
'Years elapsed',
'Sea level rise (m)', # Reorder columns
'Bruun factor (-)', dump_df = dump_df[[
'Bruun factor x SLR (m)', 'Run ID',
'Underlying trend rate (m/yr)', 'Event year',
'Underlying trend (m)', 'Years elapsed',
'Underlying + SLR (m)', 'Sea level rise (m)',
'Total movement (m)', 'Bruun factor (-)',
'Storm demand distance (m)', 'Bruun factor x SLR (m)',
'Storm demand volume (m3/m)', 'Underlying trend rate (m/yr)',
]] 'Underlying trend (m)',
'Underlying + SLR (m)',
# Sort based on maximum movement 'Total movement (m)',
dump_df = dump_df.sort_values('Total movement (m)', 'Storm demand distance (m)',
ascending=False) 'Storm demand volume (m3/m)',
]]
# Add encounter probabilities
dump_df['Encounter probability (%)'] = np.linspace( # Sort based on maximum movement
0, 100, num=n_runs + 2)[1:-1] dump_df = dump_df.sort_values('Total movement (m)',
dump_df = dump_df.set_index('Encounter probability (%)') ascending=False)
csv_name = os.path.join( # Add encounter probabilities
'diagnostics', dump_df['Encounter probability (%)'] = np.linspace(0,
'{} {} {}.csv'.format(beach_scenario, year, 100,
profile_type)) num=n_runs +
dump_df.to_csv(csv_name, float_format='%g') 2)[1:-1]
dump_df = dump_df.set_index('Encounter probability (%)')
for i, c in enumerate(dump_df.columns[3:]):
ax[i, j].plot(dump_df.index, csv_name = os.path.join(
dump_df[c],
'.',
color='#666666',
markersize=2)
ax[i, j].spines['right'].set_visible(False)
ax[i, j].spines['top'].set_visible(False)
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, 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', 'diagnostics',
f'{beach_scenario} {profile_type} scatter.png') '{} {} {}.csv'.format(beach_scenario, year, profile_type))
plt.savefig(figname, bbox_inches='tight', dpi=300) dump_df.to_csv(csv_name, float_format='%g')
plt.close(fig)
for i, c in enumerate(dump_df.columns[3:]):
ax[i, j].plot(dump_df.index,
dump_df[c],
'.',
color='#666666',
markersize=2)
ax[i, j].spines['right'].set_visible(False)
ax[i, j].spines['top'].set_visible(False)
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, 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)
figname = os.path.join(
'diagnostics', '{} {}.png'.format(beach_scenario,
profile_type))
plt.savefig(figname, bbox_inches='tight', dpi=300)
def main(): def main():

80
slr/distribution_fitting.py
Unescape Escape View File

@ -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()
Write Preview
Loading…
Cancel
Save