"Common" not failing glaciers comparison of different preprocessed gdirs options¶
In this notebook, we will compare the volume evolution, total mass change and specific mass balance of different preprocessed level 5 gdirs. For that, we always look at the largest common amount of glaciers (i.e. we compare only those glaciers that work for all preprocessed level 5 gdirs options). If you want to know about the different gdirs options that we compare, go to this section of the error analysis notebook.
If you are only interested in the final plots and a summary for preprocessing level 5 go to (after clicking scroll up to see the plot!):
- using "historical_run_output_extended" to estimate:
- regional mass change or specific mass change with qc3
- regional mass change or specific mass change: comparison qc0 vs qc3
Starting with OGGM v1.6, a small bug that occurs during the match_geod_pergla calibration (due to the use of hydro_month = 1) was corrected. This results in slightly different estimates.
In [1]:
from oggm import cfg, workflow, utils, shop
import pandas as pd
import os, glob
import numpy as np
import xarray as xr
import seaborn as sns
import matplotlib.pyplot as plt
sns.set_style("whitegrid")
cfg.initialize()
import seaborn as sns
sns.set_context('talk')
In [2]:
pd_rel_error_area_L5 = pd.read_csv('rel_error_area_statistics_for_prepro_level5_gdirs.csv', index_col=[0])
pd_working = pd.read_csv('working_rgis_for_prepro_level5_gdirs.csv', index_col='rgiid')
In [3]:
all_running_rgis = pd_working['all_running_rgis'].dropna().index.values
In [4]:
pd_geodetic = utils.get_geodetic_mb_dataframe()[utils.get_geodetic_mb_dataframe().period=='2000-01-01_2020-01-01']
# we do not look at RGI region 19 because we want to compare it to CRU!!!
pd_geodetic = pd_geodetic[pd_geodetic.reg != 19]
In [5]:
pd_geodetic_running = pd_geodetic.loc[all_running_rgis]
1. Comparison by looking at the volume evolution via "historical_run_output_extended"¶
In [6]:
rgi_reg = 'all_without_19'
pcps = ['pcp2.5', 'pcp1.6']
clims = ['CRU', 'ERA5']
dvol = {}
dvol_all = {}
print('missing combination options:')
for qc in ['qc3', 'qc0']:
for exp in ['elev_bands', 'centerlines']:
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
for pcp, clim in zip(pcps, clims):
if qc == 'qc0':
border = 80
else:
border = 160
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
try:
# todo remove this, when all preprocessed gdirs are inside of main oggm
if 'match_geod_pergla' in match:
path = '/home/www/oggm/gdirs/oggm_v1.6/'
else:
path = '/home/www/oggm/gdirs/oggm_v1.4/'
fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{pcp}/{match}/RGI62/b_{border:03d}/L5/summary/'
fs = glob.glob(fd+'historical_run_output_extended*.nc')
if rgi_reg == 'all_without_19':
for f in fs:
if 'output_extended_19' in f:
fs.remove(f)
assert len(fs) == 18
elif rgi_reg == 'all':
assert len(fs) == 19
else:
for f in fs:
if f'output_extended_{rgi_reg}' in f:
fs = [f]
break
assert len(fs) == 1
df = []
for f in fs:
ds = xr.open_dataset(f)
df.append(ds)
df = xr.concat(df, dim='rgi_id', fill_value=np.NaN)
dvol[dict_key] = df.sel(rgi_id=all_running_rgis).volume.sum(dim='rgi_id')
# do not remove those glaciers that do not work (just for comparison!!!)
dvol_all[dict_key] = df.volume.sum(dim='rgi_id')
except:
print(dict_key)
In [7]:
colors = {'no_match':'black', 'match_geod':'blue', 'match_geod_pergla': 'red', 'match_geod_pergla_massredis': 'green'}
1.1 start by only looking at qc3 to compare elev_bands vs centerlines¶
In [8]:
border = 160
qc = 'qc3'
plt.rc('font', size=18)
plt.figure(figsize=(26,14))
for exp in ['centerlines', 'elev_bands']:
for match in ['no_match', 'match_geod']: #, 'match_geod_pergla', 'match_geod_pergla_massredis']: #
for pcp, clim in zip(pcps, clims):
if exp =='centerlines':
alpha=0.5
else:
alpha = 1
if clim == 'CRU':
ls = '--'
else:
ls = '-'
try:
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
label = f'{exp}_{pcp}_{clim}_{match}'
plt.plot(dvol[dict_key].time[dvol[dict_key].time>=1900], dvol[dict_key][dvol[dict_key].time>=1900].values/1e9,
label=label,alpha=alpha,
ls=ls, color=colors[match])
except:
pass
plt.legend(title = f'{qc}_b{border}_rgi_{rgi_reg}')
ax = plt.gca()
plt.ylabel('Global volume (in km3)\n(without failing glaciers and without RGI 19)', fontsize=20)
plt.xlabel('(hydro) year', fontsize=20);
In [9]:
border = 160
qc = 'qc3'
plt.rc('font', size=18)
plt.figure(figsize=(26,14))
for exp in ['centerlines', 'elev_bands']:
for match in ['no_match', 'match_geod']: #, 'match_geod_pergla', 'match_geod_pergla_massredis']: #
for pcp, clim in zip(pcps, clims):
if exp =='centerlines':
alpha=0.5
else:
alpha = 1
if clim == 'CRU':
ls = '--'
else:
ls = '-'
try:
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
label = f'{exp}_{pcp}_{clim}_{match}'
plt.plot(dvol[dict_key].time[dvol[dict_key].time>=2000], dvol[dict_key][dvol[dict_key].time>=2000].values/1e9,
label=label, alpha=alpha,
ls=ls, color=colors[match])
except:
pass
plt.legend(title = f'{qc}_b{border}_rgi_{rgi_reg}')
ax = plt.gca()
plt.xlim([1999.5,2020.5])
plt.xticks(np.arange(2000,2021,2));
plt.ylabel('Global volume (in km3)\n(without failing glaciers and without RGI 19)', fontsize=20)
plt.xlabel('(hydro) year', fontsize=20);
plt.savefig(f'volume_change_2000_2020_for_prepro_level_gdirs_{qc}_b{border}.png', bbox_inches='tight',
pad_inches=0)
plt.savefig(f'volume_change_2000_2020_for_prepro_level_gdirs_{qc}_b{border}.pdf', bbox_inches='tight',
pad_inches=0.2)
- using centerlines instead of elev_bands results in a slightly larger volume over the 20-year time period!
1.2 now compare qc0 to qc3 and look at all match options for elev_bands only!¶
In [10]:
plt.rc('font', size=18)
plt.figure(figsize=(26,14))
alpha = 1
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
for pcp, clim in zip(pcps, clims):
for qc in ['qc0', 'qc3']:
if qc == 'qc0':
border = 80
alpha_min = 0
else:
border = 160
alpha_min = 0.5
if clim == 'CRU':
ls = '--'
else:
ls = '-'
if qc=='qc3' and 'match_geod_pergla' in match: # we don't want to look at qc3 and match_geod_pergla
pass
else:
try:
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
label = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}'
plt.plot(dvol[dict_key].time[dvol[dict_key].time>=2000], dvol[dict_key][dvol[dict_key].time>=2000].values/1e9,
label=label,alpha=alpha-alpha_min,
ls=ls, color=colors[match])
except:
print('missing: {}'.format(dict_key))
plt.legend(title = 'qc0 (no external climate quality check and correction) vs qc3')
ax = plt.gca()
plt.xlim([1999.5,2020.5])
plt.xticks(np.arange(2000,2021,2));
plt.ylabel('Global volume (in km3)\n(without failing glaciers and without RGI 19)', fontsize=20)
plt.xlabel('(hydro) year', fontsize=20);
plt.savefig(f'volume_change_2000_2020_for_prepro_level_gdirs_qc0_qc3.png', bbox_inches='tight',
pad_inches=0)
plt.savefig(f'volume_change_2000_2020_for_prepro_level_gdirs_qc0_qc3.pdf', bbox_inches='tight',
pad_inches=0.2)
when looking at all glaciers that work for all level 5 preprocessed gdirs with elev_bands:
- in 2020:
- no_match results in smallest world-wide volume, match_geod in largest and match_geod_pergla has a slightly smaller volume than match_geod
- the largest differences in the volume between CRU and ERA5 occur between 2007 and 2017 for all preprocessed gdirs (for some they diverge earlier for others they come together later ...)
in 2000:
- no_match has the largest volume, followed by match_geod and match_geod_pergla has the smallest volume
match_geod_pergla_massredis is equal to match_geod_pergla in 2000, but over time, using a mass redistribution curve instead of the SIA results in very slightly higher volume estimates
strongest difference between qc0 and qc3: (only looking at CRU & elev_bands here)
- in the absolute volume 2020: for match_geod (although in terms of volume change they are almost identical, see dmdt plot below)
- in the relative volume change: for no_match (qc0 has a much more negative volume change than qc3, see dmdt plot below!)
If you are interested in how these method match options produce different volume projections in the future, you can check out this notebook (only for CRU comparing in total 3 options)!
In [ ]:
1.3 compare volume changes of common running glaciers to all running glaciers:¶
- a) using only common running glaciers that work for all preprocessed gdirs (as before) to
- b) using all glaciers that work for each preprocessed gdirs. This means there are different amount of glaciers per option!
In [11]:
border = 80 #160
qc = 'qc0'
plt.rc('font', size=18)
plt.figure(figsize=(26,14))
for exp in ['centerlines', 'elev_bands']:
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
for pcp, clim in zip(pcps, clims):
if exp =='centerlines':
alpha=0.5
else:
alpha = 1
if clim == 'CRU':
ls = '--'
else:
ls = '-'
try:
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
label = f'{exp}_{pcp}_{clim}_{match}_all'
plt.plot(dvol_all[dict_key].time[dvol_all[dict_key].time>=2000], dvol_all[dict_key][dvol_all[dict_key].time>=2000].values/1e9,
label=label,alpha=alpha,
ls = ls, color = colors[match], linewidth=4)
label = f'{exp}_{pcp}_{clim}_{match}_only_common_running_glaciers'
plt.plot(dvol[dict_key].time[dvol[dict_key].time>=2000], dvol[dict_key][dvol[dict_key].time>=2000].values/1e9,
label=label,alpha=alpha,
ls = ls, color = colors[match], linewidth=1)
except:
print(dict_key)
plt.legend( title = f'{qc}_b{border}_rgi_{rgi_reg}', fontsize=16)
ax = plt.gca()
plt.xlim([1999.5,2020.5])
plt.xticks(np.arange(2000,2021,2));
plt.ylabel('Global volume (in km3)\n(without failing glaciers and without RGI 19)', fontsize=20)
plt.xlabel('(hydro) year', fontsize=20);