Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add uncertainty budget plot #8

Closed
wants to merge 1 commit into from
Closed

Add uncertainty budget plot #8

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 51 additions & 13 deletions step4_plot_regional_sealevel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@

from config import settings
from tide_gauge_locations import extract_site_info, find_nearest_station_id
from slr_pkg import abbreviate_location_name # found in __init.py__
from slr_pkg import abbreviate_location_name, choose_montecarlo_dir # found in __init.py__
from surge import tide_gauge_library as tgl
from directories import read_dir, makefolder
from plotting_libraries import location_string, scenario_string, \
Expand All @@ -29,12 +29,13 @@ def compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl):
:return: years, scenario uncertainty, model uncertainty, internal
variability
"""
allpmid = np.zeros([3, 94])
allunc = np.zeros([3, 94])
nyrs = np.arange(2007, settings["projection_end_year"] + 1).size
allpmid = np.zeros([3, nyrs])
allunc = np.zeros([3, nyrs])

# Estimate internal variability from de-trended gauge data
tg_years_arr = np.array(tg_years, dtype='int')
vals = np.ma.masked_values(tg_amsl, -99999.)
vals = np.ma.masked_invalid(tg_amsl)
IntV = compute_variability(tg_years_arr / 1000., vals / 1.)

# Get the values of the multi-index: years and percentile values
Expand All @@ -53,7 +54,7 @@ def compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl):
# Model Uncertainty (already expressed as 90% confidence interval)
UncM = np.mean(allunc, axis=0)

return years, UncS, UncM, IntV
return years, np.squeeze(UncS), np.squeeze(UncM), IntV


def compute_variability(x, y, factor=1.645):
Expand All @@ -64,10 +65,8 @@ def compute_variability(x, y, factor=1.645):
:param factor: multiplication factor
:return: internal variability component of uncertainty
"""
mask = np.ma.getmask(y)
index = np.where(mask is True)
new_x = np.delete(x, index)
new_y = np.delete(y, index)
new_x = x[~y.mask]
new_y = y[~y.mask]
fit = np.polyfit(new_x, new_y, 1)
fit_data = new_x * fit[0] + fit[1]
stdev = np.std(new_y - fit_data)
Expand Down Expand Up @@ -711,6 +710,37 @@ def plot_figure_seven(g_df_list, r_df_list, site_name, scenarios, fig_dir):
plt.close()


def plot_figure_eight(df_r_list, tg_years, tg_amsl, site_name, scenarios, fig_dir):
years, UncS, UncM, IntV = compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl)

IntV = np.full(UncS.shape, IntV)

UncT = np.sum([UncS, UncM, IntV], axis=0)

fig = plt.figure(figsize=(4.5, 4.5))

ax = fig.add_subplot(111)
matplotlib.rcParams['font.size'] = 9

ax.fill_between(years, 0, UncM/UncT, color='blue', label='Model')
ax.fill_between(years, UncM/UncT, (UncM+UncS)/UncT, color='green', label='Scenario')
ax.fill_between(years, (UncM + UncS)/UncT, 1, color='orange', label='Variability')

ax.legend(loc='lower right', fancybox=True, framealpha=1)
ax.set_title(f'{location_string(site_name)}')
ax.set_xlabel('Year')
ax.set_ylabel('Fraction of total variance')

ax.set_xlim(2000, years[-1])
ax.set_ylim(0, 1)

fig.tight_layout()
outfile = f'{fig_dir}08_{site_name}_uncertainty_budget.png'
plt.savefig(outfile, dpi=200, format='png')
plt.show()
plt.close()


def plot_tg_data(ax, nflag, flag, tg_years, non_missing, tg_amsl, tg_name):
"""
Plot the annual mean sea levels from the tide gauge data.
Expand Down Expand Up @@ -787,19 +817,21 @@ def read_IPCC_AR5_Levermann_proj(scenarios, refname='sum'):
print('running function read_IPCC_AR5_Levermann_proj')

# Directory of Monte Carlo time series for new projections
mcdir = settings["montecarlodir"]
mcdir = choose_montecarlo_dir()

nyrs = np.arange(2007, settings["projection_end_year"] + 1).size

ar5_low = []
ar5_mid = []
ar5_upp = []

for sce in scenarios:
reflow = iris.load_cube(
os.path.join(mcdir, sce + '_' + refname + 'lower.nc')).data
os.path.join(mcdir, sce + '_' + refname + 'lower.nc'))[:nyrs].data
refmid = iris.load_cube(
os.path.join(mcdir, sce + '_' + refname + 'mid.nc')).data
os.path.join(mcdir, sce + '_' + refname + 'mid.nc'))[:nyrs].data
refupp = iris.load_cube(
os.path.join(mcdir, sce + '_' + refname + 'upper.nc')).data
os.path.join(mcdir, sce + '_' + refname + 'upper.nc'))[:nyrs].data

ar5_low.append(reflow)
ar5_mid.append(refmid)
Expand Down Expand Up @@ -964,6 +996,12 @@ def main():
# Same style as Figure 4 Howard and Palmer (2019)
plot_figure_seven(g_df_list, r_df_list, df_loc, rcp_scenarios,
sealev_fdir)

# Figure 8
# Uncertainty budget for local sea level projections,
# based on Hawkins & Sutton (2009)
plot_figure_eight(r_df_list, tg_years, tg_amsl, df_loc,
rcp_scenarios, sealev_fdir)


if __name__ == '__main__':
Expand Down