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!):
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.
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')
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')
all_running_rgis = pd_working['all_running_rgis'].dropna().index.values
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]
pd_geodetic_running = pd_geodetic.loc[all_running_rgis]
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)
colors = {'no_match':'black', 'match_geod':'blue', 'match_geod_pergla': 'red', 'match_geod_pergla_massredis': 'green'}
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);
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)
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 2000:
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)
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)!
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);