"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!):

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)

  • 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)
missing: elev_bands_pcp1.6_ERA5_no_match_qc0_b80_rgi_all_without_19
missing: elev_bands_pcp1.6_ERA5_match_geod_qc0_b80_rgi_all_without_19

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);
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
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