From e86b7da54610706be06cf7faecc6043e21b459c5 Mon Sep 17 00:00:00 2001 From: mo-gregmunday Date: Wed, 12 Jun 2024 14:13:59 +0100 Subject: [PATCH 1/3] Fixed uncertainty budget calculations and added plot --- step4_plot_regional_sealevel.py | 64 ++++++++++++++++++++++++++------- 1 file changed, 51 insertions(+), 13 deletions(-) diff --git a/step4_plot_regional_sealevel.py b/step4_plot_regional_sealevel.py index 55b34e9..2507a41 100644 --- a/step4_plot_regional_sealevel.py +++ b/step4_plot_regional_sealevel.py @@ -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, \ @@ -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 @@ -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): @@ -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) @@ -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. @@ -787,7 +817,9 @@ 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 = [] @@ -795,11 +827,11 @@ def read_IPCC_AR5_Levermann_proj(scenarios, refname='sum'): 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) @@ -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__': From 269d906c41506d9851bb1d6fdef9bdd339b2e688 Mon Sep 17 00:00:00 2001 From: mo-gregmunday Date: Wed, 12 Jun 2024 14:24:20 +0100 Subject: [PATCH 2/3] Update with old catch --- step4_plot_regional_sealevel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/step4_plot_regional_sealevel.py b/step4_plot_regional_sealevel.py index 2507a41..8034138 100644 --- a/step4_plot_regional_sealevel.py +++ b/step4_plot_regional_sealevel.py @@ -36,6 +36,7 @@ def compute_uncertainties(df_r_list, scenarios, tg_years, tg_amsl): # Estimate internal variability from de-trended gauge data tg_years_arr = np.array(tg_years, dtype='int') vals = np.ma.masked_invalid(tg_amsl) + vals = np.ma.masked_values(vals, -99999.) IntV = compute_variability(tg_years_arr / 1000., vals / 1.) # Get the values of the multi-index: years and percentile values @@ -737,7 +738,6 @@ def plot_figure_eight(df_r_list, tg_years, tg_amsl, site_name, scenarios, fig_di fig.tight_layout() outfile = f'{fig_dir}08_{site_name}_uncertainty_budget.png' plt.savefig(outfile, dpi=200, format='png') - plt.show() plt.close() From ba4061f4d5c63cbe70992a441fcdfe429bdf13a8 Mon Sep 17 00:00:00 2001 From: mo-gregmunday Date: Thu, 13 Jun 2024 10:33:08 +0100 Subject: [PATCH 3/3] Tide gauge as lines --- step4_plot_regional_sealevel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/step4_plot_regional_sealevel.py b/step4_plot_regional_sealevel.py index c412932..7e4ab47 100644 --- a/step4_plot_regional_sealevel.py +++ b/step4_plot_regional_sealevel.py @@ -758,11 +758,11 @@ def plot_tg_data(ax, nflag, flag, tg_years, non_missing, tg_amsl, tg_name): print(f'Tide gauge data has been flagged for attention - ' + f'{tg_years[(flag & non_missing)]}') ax.plot(tg_years[(flag & non_missing)], tg_amsl[(flag & non_missing)], - marker='o', mec='black', mfc='None', - markersize=3, linestyle='None', label='TG flagged') + lw=1, color='black', mfc='None', + linestyle='None', label='TG flagged') if nflag < len(flag): ax.plot(tg_years[(~flag & non_missing)], - tg_amsl[(~flag & non_missing)], 'ko', markersize=3, + tg_amsl[(~flag & non_missing)], lw=1, color='black', label=f'{location_string(tg_name)} TG')