"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')
2022-07-25 08:45:28: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file. 2022-07-25 08:45:28: oggm.cfg: Multiprocessing switched OFF according to the parameter file. 2022-07-25 08:45:28: oggm.cfg: Multiprocessing: using all available processors (N=32) 2022-07-25 08:45:28: oggm.utils: Checking the download verification file checksum...
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')
/home/users/lschuster/.local/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3444: DtypeWarning: Columns (1,2,3,4) have mixed types.Specify dtype option on import or set low_memory=False. exec(code_obj, self.user_global_ns, self.user_ns)
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]
2022-07-25 08:45:30: oggm.utils: No known hash for cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb/hugonnet_2021_ds_rgi60_pergla_rates_10_20_worldwide_filled.hdf
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)
missing combination options: elev_bands_pcp2.5_CRU_match_geod_pergla_qc3_b160_rgi_all_without_19 elev_bands_pcp1.6_ERA5_match_geod_pergla_qc3_b160_rgi_all_without_19 elev_bands_pcp2.5_CRU_match_geod_pergla_massredis_qc3_b160_rgi_all_without_19 elev_bands_pcp1.6_ERA5_match_geod_pergla_massredis_qc3_b160_rgi_all_without_19 centerlines_pcp2.5_CRU_match_geod_qc3_b160_rgi_all_without_19 centerlines_pcp1.6_ERA5_match_geod_qc3_b160_rgi_all_without_19 centerlines_pcp2.5_CRU_match_geod_pergla_qc3_b160_rgi_all_without_19 centerlines_pcp1.6_ERA5_match_geod_pergla_qc3_b160_rgi_all_without_19 centerlines_pcp2.5_CRU_match_geod_pergla_massredis_qc3_b160_rgi_all_without_19 centerlines_pcp1.6_ERA5_match_geod_pergla_massredis_qc3_b160_rgi_all_without_19 elev_bands_pcp1.6_ERA5_no_match_qc0_b80_rgi_all_without_19 elev_bands_pcp1.6_ERA5_match_geod_qc0_b80_rgi_all_without_19 centerlines_pcp2.5_CRU_no_match_qc0_b80_rgi_all_without_19 centerlines_pcp1.6_ERA5_no_match_qc0_b80_rgi_all_without_19 centerlines_pcp2.5_CRU_match_geod_qc0_b80_rgi_all_without_19 centerlines_pcp1.6_ERA5_match_geod_qc0_b80_rgi_all_without_19 centerlines_pcp2.5_CRU_match_geod_pergla_qc0_b80_rgi_all_without_19 centerlines_pcp1.6_ERA5_match_geod_pergla_qc0_b80_rgi_all_without_19 centerlines_pcp2.5_CRU_match_geod_pergla_massredis_qc0_b80_rgi_all_without_19 centerlines_pcp1.6_ERA5_match_geod_pergla_massredis_qc0_b80_rgi_all_without_19
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)