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')
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...
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)
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]
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
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)
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
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)
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:
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)!
to
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
# first get W5E5 volume changes (only available for qc0, elev_band, match_geod_pergla, border=160)
pcps_w5e5 = ['2.6', 'via_winterprcp']
clims_w5e5 = ['GSWP3_W5E5', 'GSWP3_W5E5']
for qc in ['qc0']:
for exp in ['elev_bands']:
for match in ['match_geod_pergla']: #
for pcp, clim in zip(pcps_w5e5, clims_w5e5):
#if qc == 'qc0':
# border = 80
#else:
border = 160
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
#try:
path = '/home/www/fmaussion/runs/new_gdirs/oggm_v1.6'
#fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{pcp}/{match}/RGI62/b_{border:03d}/L5/summary/'
if pcp == '2.6':
matchi = match + '_2.6'
else:
matchi = match + '_winterprcp'
fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{matchi}/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)
clims_all = ['CRU', 'ERA5', 'GSWP3_W5E5', 'GSWP3_W5E5']
pcps_all = ['pcp2.5', 'pcp1.6','2.6', 'via_winterprcp']
plt.rc('font', size=18)
plt.figure(figsize=(26,14))
alpha = 1
for clim, pcp in zip(clims_all, pcps_all):
for qc in ['qc0']: # , 'qc3']:
for match in ['match_geod_pergla']: #, 'match_geod_pergla_massredis']: #
if 'W5E5' in clim:
border = 160
color = 'blue'
else:
border= 80
#alpha_min = 0.7
if clim == 'CRU':
color = 'red'
elif clim == 'ERA5':
color = 'green'
if 'winterprcp' in pcp:
ls = '--'
else:
ls = '-'
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=color)
except:
print('missing: {}'.format(dict_key))
plt.legend(title = 'Climate dataset & prcp. fac choice comparison')
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_climate_pcp_comparison.png', bbox_inches='tight',
pad_inches=0)
plt.savefig(f'volume_change_2000_2020_for_prepro_level_gdirs_climate_pcp_comparison.pdf', bbox_inches='tight',
pad_inches=0.2)
period = '2000-01-01_2020-01-01'
#dfz = pd.read_csv(utils.get_demo_file('zemp_ref_2006_2016.csv'), index_col=0)
#dfh['dmdt_zemp'] = dfz.SMB.values * 1000
#dfh['dmdt_zemp_err'] = dfz.SMB_err.values * 1000
dfh = pd.read_csv(utils.get_demo_file('table_hugonnet_regions_10yr_20yr_ar6period.csv'), index_col=0)
dfh = dfh.loc[dfh.period == period]
dfh.index = ['{:02d}'.format(int(rgi_reg)) for rgi_reg in dfh.index]
dfh = dfh.drop('19')
run = False
pcps = ['pcp2.5', 'pcp1.6']
clims = ['CRU', 'ERA5']
if run:
for qc in ['qc0', 'qc3']:
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_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}'
# todo remove this, when all preprocessed gdirs are inside of main oggm
if 'match_geod_pergla' in match and qc == 'qc0':
path = '/home/www/oggm/gdirs/oggm_v1.6/'
elif 'match_geod_pergla' in match and qc == 'qc3':
break
elif 'match_geod_pergla' in match and exp == 'centerlines':
break
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/'
for rgi_reg in np.arange(1, 19):
# for comparison, we do not want to look at Antarctica (RGI region19)
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
rgi_reg_int = rgi_reg.copy()
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
rgi_reg = '{:02d}'.format(rgi_reg)
try:
df = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
dfs = pd.read_csv(fd + f'glacier_statistics_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
except FileNotFoundError:
#print('Not here:',rgi_reg)
continue
df = df.dropna(axis=0, how='all')
# just choose those glaciers that work for all prepro gdirs types!
df = df[all_running_rgis_reg]
# check if there are no np.NaNs (otherwise sth. is wrong with all_running_rgis_reg)
assert ~np.any(df.isna())
#df = df.dropna(axis=1, how='all')
# odf = pd.DataFrame(df.loc[2006:2016].mean(), columns=['SMB'])
# odf['AREA'] = dfs.rgi_area_km2
# dfh.loc[rgi_reg, 'AREA_OGGM'] = odf['AREA'].sum()
# dfh.loc[rgi_reg, 'SMB_OGGM'] = np.average(odf['SMB'], weights=odf['AREA']) / 1000
odf = pd.DataFrame(df.loc[2000:].mean(), columns=['SMB'])
odf['AREA'] = dfs.rgi_area_km2
dfh.loc[rgi_reg, f'dmdtda_OGGM_{dict_key_short}'] = np.average(odf['SMB'], weights=odf['AREA'])
#odir_5 = odir.replace('L3/', 'L5/')
with xr.open_dataset(fd + f'historical_run_output_extended_{rgi_reg}.nc') as ds:
ds = ds[['volume', 'area']].sum(dim='rgi_id')
vol_ts = ds.volume.to_series()
area_ts = ds.area.to_series()
# dmdt is in kg per year *10e-12
dfh.loc[rgi_reg, f'dmdt_OGGM_{dict_key_short}'] = (vol_ts.loc[2019] - vol_ts.loc[2000]) * cfg.PARAMS['ice_density'] * 1e-12 / 20
dfh.loc[rgi_reg, f'dmdtda_dyna_OGGM_{dict_key_short}'] = (vol_ts.loc[2019] - vol_ts.loc[2000]) / area_ts.loc[2000] * cfg.PARAMS['ice_density'] / 20
# also get W5E5 volume changes (only available for qc0, elev_band, match_geod_pergla, border=160)
pcps_w5e5 = ['2.6', 'via_winterprcp']
clims_w5e5 = ['GSWP3_W5E5', 'GSWP3_W5E5']
for qc in ['qc0']:
for exp in ['elev_bands']:
for match in ['match_geod_pergla']: #
for pcp, clim in zip(pcps_w5e5, clims_w5e5):
#if qc == 'qc0':
# border = 80
#else:
border = 160
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}'
#try:
path = '/home/www/fmaussion/runs/new_gdirs/oggm_v1.6'
#fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{pcp}/{match}/RGI62/b_{border:03d}/L5/summary/'
if pcp == '2.6':
matchi = match + '_2.6'
else:
matchi = match + '_winterprcp'
fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{matchi}/RGI62/b_{border:03d}/L5/summary/'
for rgi_reg in np.arange(1, 19):
# for comparison, we do not want to look at Antarctica (RGI region19)
dict_key = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
rgi_reg_int = rgi_reg.copy()
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
rgi_reg = '{:02d}'.format(rgi_reg)
try:
df = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
dfs = pd.read_csv(fd + f'glacier_statistics_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
except FileNotFoundError:
#print('Not here:',rgi_reg)
continue
df = df.dropna(axis=0, how='all')
# just choose those glaciers that work for all prepro gdirs types!
df = df[all_running_rgis_reg]
# check if there are no np.NaNs (otherwise sth. is wrong with all_running_rgis_reg)
assert ~np.any(df.isna())
#df = df.dropna(axis=1, how='all')
# odf = pd.DataFrame(df.loc[2006:2016].mean(), columns=['SMB'])
# odf['AREA'] = dfs.rgi_area_km2
# dfh.loc[rgi_reg, 'AREA_OGGM'] = odf['AREA'].sum()
# dfh.loc[rgi_reg, 'SMB_OGGM'] = np.average(odf['SMB'], weights=odf['AREA']) / 1000
odf = pd.DataFrame(df.loc[2000:].mean(), columns=['SMB'])
odf['AREA'] = dfs.rgi_area_km2
dfh.loc[rgi_reg, f'dmdtda_OGGM_{dict_key_short}'] = np.average(odf['SMB'], weights=odf['AREA'])
#odir_5 = odir.replace('L3/', 'L5/')
with xr.open_dataset(fd + f'historical_run_output_extended_{rgi_reg}.nc') as ds:
ds = ds[['volume', 'area']].sum(dim='rgi_id')
vol_ts = ds.volume.to_series()
area_ts = ds.area.to_series()
# dmdt is in kg per year *10e-12
dfh.loc[rgi_reg, f'dmdt_OGGM_{dict_key_short}'] = (vol_ts.loc[2019] - vol_ts.loc[2000]) * cfg.PARAMS['ice_density'] * 1e-12 / 20
dfh.loc[rgi_reg, f'dmdtda_dyna_OGGM_{dict_key_short}'] = (vol_ts.loc[2019] - vol_ts.loc[2000]) / area_ts.loc[2000] * cfg.PARAMS['ice_density'] / 20
dfh.to_csv(f'dmdtda_dmdt_for_prepro_level_5_gdirs_new.csv')
else:
# new with (partly) preprocessed gdirs
dfh = pd.read_csv(f'dmdtda_dmdt_for_prepro_level_5_gdirs_new.csv', index_col=[0])
# old with wrong preprocessed gdirs: dfh = pd.read_csv(f'dmdtda_dmdt_for_prepro_level_5_gdirs.csv', index_col=[0])
dfh.index = ['{:02d}'.format(int(rgi_reg)) for rgi_reg in dfh.index]
dfh.columns
Index(['period', 'dmdt', 'err_dmdt', 'dmdtda', 'err_dmdtda', 'dmdtda_OGGM_elev_bands_pcp2.5_CRU_no_match_qc0_b80', 'dmdt_OGGM_elev_bands_pcp2.5_CRU_no_match_qc0_b80', 'dmdtda_dyna_OGGM_elev_bands_pcp2.5_CRU_no_match_qc0_b80', 'dmdtda_OGGM_elev_bands_pcp2.5_CRU_match_geod_qc0_b80', 'dmdt_OGGM_elev_bands_pcp2.5_CRU_match_geod_qc0_b80', 'dmdtda_dyna_OGGM_elev_bands_pcp2.5_CRU_match_geod_qc0_b80', 'dmdtda_OGGM_elev_bands_pcp2.5_CRU_match_geod_pergla_qc0_b80', 'dmdt_OGGM_elev_bands_pcp2.5_CRU_match_geod_pergla_qc0_b80', 'dmdtda_dyna_OGGM_elev_bands_pcp2.5_CRU_match_geod_pergla_qc0_b80', 'dmdtda_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_qc0_b80', 'dmdt_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_qc0_b80', 'dmdtda_dyna_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_qc0_b80', 'dmdtda_OGGM_elev_bands_pcp2.5_CRU_match_geod_pergla_massredis_qc0_b80', 'dmdt_OGGM_elev_bands_pcp2.5_CRU_match_geod_pergla_massredis_qc0_b80', 'dmdtda_dyna_OGGM_elev_bands_pcp2.5_CRU_match_geod_pergla_massredis_qc0_b80', 'dmdtda_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_massredis_qc0_b80', 'dmdt_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_massredis_qc0_b80', 'dmdtda_dyna_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_massredis_qc0_b80', 'dmdtda_OGGM_elev_bands_pcp2.5_CRU_no_match_qc3_b160', 'dmdt_OGGM_elev_bands_pcp2.5_CRU_no_match_qc3_b160', 'dmdtda_dyna_OGGM_elev_bands_pcp2.5_CRU_no_match_qc3_b160', 'dmdtda_OGGM_elev_bands_pcp1.6_ERA5_no_match_qc3_b160', 'dmdt_OGGM_elev_bands_pcp1.6_ERA5_no_match_qc3_b160', 'dmdtda_dyna_OGGM_elev_bands_pcp1.6_ERA5_no_match_qc3_b160', 'dmdtda_OGGM_elev_bands_pcp2.5_CRU_match_geod_qc3_b160', 'dmdt_OGGM_elev_bands_pcp2.5_CRU_match_geod_qc3_b160', 'dmdtda_dyna_OGGM_elev_bands_pcp2.5_CRU_match_geod_qc3_b160', 'dmdtda_OGGM_elev_bands_pcp1.6_ERA5_match_geod_qc3_b160', 'dmdt_OGGM_elev_bands_pcp1.6_ERA5_match_geod_qc3_b160', 'dmdtda_dyna_OGGM_elev_bands_pcp1.6_ERA5_match_geod_qc3_b160', 'dmdtda_OGGM_centerlines_pcp2.5_CRU_no_match_qc3_b160', 'dmdt_OGGM_centerlines_pcp2.5_CRU_no_match_qc3_b160', 'dmdtda_dyna_OGGM_centerlines_pcp2.5_CRU_no_match_qc3_b160', 'dmdtda_OGGM_centerlines_pcp1.6_ERA5_no_match_qc3_b160', 'dmdt_OGGM_centerlines_pcp1.6_ERA5_no_match_qc3_b160', 'dmdtda_dyna_OGGM_centerlines_pcp1.6_ERA5_no_match_qc3_b160', 'dmdtda_OGGM_elev_bands_2.6_GSWP3_W5E5_match_geod_pergla_qc0_b160', 'dmdt_OGGM_elev_bands_2.6_GSWP3_W5E5_match_geod_pergla_qc0_b160', 'dmdtda_dyna_OGGM_elev_bands_2.6_GSWP3_W5E5_match_geod_pergla_qc0_b160', 'dmdtda_OGGM_elev_bands_via_winterprcp_GSWP3_W5E5_match_geod_pergla_qc0_b160', 'dmdt_OGGM_elev_bands_via_winterprcp_GSWP3_W5E5_match_geod_pergla_qc0_b160', 'dmdtda_dyna_OGGM_elev_bands_via_winterprcp_GSWP3_W5E5_match_geod_pergla_qc0_b160'], dtype='object')
# get those geodetic estimates of the common running glaciers and that do the mean!
dmdtda_working_glaciers_geods = {}
dmdt_working_glaciers_geods = {}
for reg in np.arange(1,19,1):
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == reg]
dmdtda_working_glaciers_geod = np.average(pd_geodetic_running_reg.dmdtda*1e3, weights=pd_geodetic_running_reg.area)
dmdt_working_glaciers_geod = (pd_geodetic_running_reg.dmdtda*1e3*pd_geodetic_running_reg.area).sum()*1e-12
dmdtda_working_glaciers_geods[reg] = dmdtda_working_glaciers_geod
dmdt_working_glaciers_geods[reg] = dmdt_working_glaciers_geod
pd_working_glaciers_geod =pd.DataFrame([dmdt_working_glaciers_geods,dmdtda_working_glaciers_geods],
index=['dmdt_geodetic_only_running_glaciers', 'dmdtda_geodetic_only_running_glaciers']).astype(float)
dfh['dmdtda_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_qc0_b80']
01 -775.111072 02 -521.323876 03 -291.469884 04 -648.097569 05 -411.157894 06 -823.917118 07 -296.225023 08 -565.072236 09 -202.368134 10 -524.435265 11 -805.371476 12 -508.414598 13 -234.597791 14 -136.422764 15 -472.069305 16 -390.949152 17 -734.417725 18 -556.277795 Name: dmdtda_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_qc0_b80, dtype: float64
border = 160
qc = 'qc3'
plt.figure(figsize=(26,26))
n=17
dfh = pd.read_csv('dmdtda_dmdt_for_prepro_level_5_gdirs_new.csv', index_col=[0])
dfh.index = ['{:02d}'.format(int(rgi_reg)) for rgi_reg in dfh.index]
dfh['dmdt_geodetic_only_running_glaciers'] = pd_working_glaciers_geod.T['dmdt_geodetic_only_running_glaciers'].values
dfh['dmdtda_geodetic_only_running_glaciers'] = pd_working_glaciers_geod.T['dmdtda_geodetic_only_running_glaciers'].values
###
###
ax = plt.subplot2grid((n, n), (0, 0), colspan=8, rowspan=8)
#plt.figure(figsize=(20,10))
ax = plt.gca()
#f, ax = plt.subplots()
alphas = {'elev_bands': 1, 'centerlines':0.5}
markers = {'ERA5': 'o', 'CRU': 'x'}
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):
try:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdt_OGGM_{dict_key_short}', marker=markers[clim],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp],
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
dfh.plot(ax=ax, y='dmdt', yerr='err_dmdt', marker='s', linestyle='none', markersize=12,
label = 'geodetic_observation_Hugonnet2021', color = 'orange', alpha = 0.5)
dfh.plot(ax=ax, y='dmdt_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation only common running glaciers', color = 'violet', alpha = 1, markeredgewidth=5)
plt.ylabel(r'regional total mass change (dmdt, Gt year$^{-1}$)')
plt.xlabel('RGI region')
plt.xticks(np.arange(0,18,1), dfh.index.values)
f = ax.get_legend()
f.remove()
###
ax = plt.subplot2grid((n, n), (9, 0), colspan=8, rowspan=8)
#f, ax = plt.subplots()
alphas = {'elev_bands': 1, 'centerlines':0.5}
markers = {'ERA5': 'o', 'CRU': 'x'}
dfh.plot(ax=ax, y='dmdtda', yerr='err_dmdtda', marker='s', linestyle='none', markersize=12, linewidth=4,
label = 'geodetic observation (Hugonnet et al., 2021)', color = 'orange', alpha = 0.7)
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):
try:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdtda_OGGM_{dict_key_short}', marker=markers[clim],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp],
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
#t=ax.get_legend_handles_labels()
#t[0] = [t[0][-1]].append(t[:-1])
dfh.plot(ax=ax, y='dmdtda_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation only common running glaciers', color = 'violet', alpha = 1, markeredgewidth=5)
plt.ylabel(r'area-weighted specific mass balance (dmdtda, kg m$^{-2}$ year$^{-1}$)')
plt.xlabel('RGI region')
plt.xticks(np.arange(0, 18, 1), dfh.index.values)
f = ax.get_legend()
f.remove()
plt.title('area-weighted mean MB from "fixed_geometry_mass_balance"')
##
#plt.figure(figsize=(20,10))
ax = plt.subplot2grid((n, n), (9, 9), colspan=8, rowspan=8)
#f, ax = plt.subplots()
alphas = {'elev_bands': 1, 'centerlines':0.5}
markers = {'ERA5': 'o', 'CRU': 'x'}
bu=[]
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):
try:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdtda_dyna_OGGM_{dict_key_short}', marker=markers[clim],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp],
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
dfh.plot(ax=ax, y='dmdtda', yerr='err_dmdtda', marker='s', linestyle='none', markersize=12,
label = 'geodetic observation (Hugonnet et al., 2021)', color = 'orange', alpha = 0.7)
dfh.plot(ax=ax, y='dmdtda_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation (only common running glaciers)', color = 'violet', alpha = 1, markeredgewidth=5)
plt.ylabel(r'dynamic specific mass balance (dmdtda, kg m$^{-2}$ year$^{-1}$)')
plt.xlabel('RGI region')
plt.title('using volume changes from "historical_run_output_extended"')
plt.xticks(np.arange(0,18,1), dfh.index.values);
plt.ylim([-1600,-50])
ax.legend(framealpha=0.5, ncol=1, loc=3, bbox_to_anchor=(0.17,1.3),
title=f'mean of 2000-2020: {qc}, border {border} (only glaciers that are\nrunning for all preprocessed gdirs used, i.e. {len(all_running_rgis)*100/len(pd_working):0.1f}%)')
#ax.legend()
######
#plt.tight_layout()
ax = plt.subplot2grid((n, n), (1, 9), rowspan=6, colspan=1)
rgi_index = ['{:02d}'.format(int(rgi_reg)) for rgi_reg in np.arange(1,19)]
dfh.loc['all_without_19', 'dmdt'] = dfh.loc[rgi_index].dmdt.sum()
dfh.loc['all_without_19', 'err_dmdt'] = dfh.loc[rgi_index].err_dmdt.sum()
dfh.loc['all_without_19', 'period'] = dfh.period[0]
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):
try:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}'
dfh.loc['all_without_19', f'dmdt_OGGM_{dict_key_short}'] = dfh.loc[rgi_index, f'dmdt_OGGM_{dict_key_short}'].sum()
except:
pass
###
#plt.figure(figsize=(3,6))
ax = plt.gca()
#ax.errorbar(['geodetic observation (Hugonnet et al., 2021)'],
# df_dmdt_all.dmdt, yerr=df_dmdt_all.err_dmdt,
# marker='s', color='orange', alpha=0.8)
df_dmdt_all = dfh.loc['all_without_19'][1:].dropna()
plt.axhline(df_dmdt_all.dmdt, color='orange', alpha = 0.3)
plt.axhspan(df_dmdt_all.dmdt-df_dmdt_all.err_dmdt,
df_dmdt_all.dmdt + df_dmdt_all.err_dmdt, alpha = 0.1, color='orange')
#dmdtda_working_glaciers_geods['all_without_19'] = np.average(pd_geodetic_running.dmdtda, weights=pd_geodetic_running.area)*1e3
plt.axhline((pd_geodetic_running.dmdtda*1e3*pd_geodetic_running.area).sum()*1e-12,
color='violet', alpha=0.5, linewidth=5)
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):
try:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
ax.errorbar(x=[f'{exp}_{clim}_{pcp}_{match}'], # [0]
y=df_dmdt_all[f'dmdt_OGGM_{dict_key_short}'],
marker=markers[clim],
linestyle='none', markersize=8, color=colors[match], alpha=alphas[exp],
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
plt.ylabel(r'world-wide total mass-change (dmdt, Gt year$^{-1}$)')
plt.xticks(ticks=[3.5],labels=['all without RGI 19'])
plt.xlim([-1,9.5])
plt.grid(axis='x')
plt.savefig(f'dmdtda_dmdt_for_prepro_level_gdirs_{qc}_b{border}.png', bbox_inches='tight', pad_inches=0)
plt.savefig(f'dmdtda_dmdt_for_prepro_level_gdirs_{qc}_b{border}.pdf', bbox_inches='tight', pad_inches=0.2)
normally match_geod should perfectly match the geodetic estimates, however, here we compare the "common" running glaciers of all preprocessed gdirs. Hence, match_geod does not always perfectly match even with the fixed geometry MB estimates
the dynamic SMB matches less good (specifically visible for match_geod) because the changing area was neglected during the MB calibration.
Somehow naturally, the no_match option is most different to the geodetic estimates (as here only the direct glaciological reference data are used):
For the geodetic observation, we distinguish between geodetic observations for all RGI glaciers of that region or geodetic estimates for only the common running glaciers
dfh = pd.read_csv('dmdtda_dmdt_for_prepro_level_5_gdirs_new.csv', index_col=[0])
dfh.index = ['{:02d}'.format(int(rgi_reg)) for rgi_reg in dfh.index]
dfh['dmdt_geodetic_only_running_glaciers'] = pd_working_glaciers_geod.T['dmdt_geodetic_only_running_glaciers'].values
dfh['dmdtda_geodetic_only_running_glaciers'] = pd_working_glaciers_geod.T['dmdtda_geodetic_only_running_glaciers'].values
markers = {'ERA5': 'o', 'CRU': 'x', 'GSWP3_W5E5_cte':'v', 'GSWP3_W5E5_var':'^'}
# only CRU and 'elev_bands' available for qc0 at the moment!!!
exp = 'elev_bands'
plt.figure(figsize=(26,26))
n=17
ax = plt.subplot2grid((n, n), (0, 0), colspan=8, rowspan=8)
ax = plt.gca()
clims_all = ['CRU', 'ERA5', 'GSWP3_W5E5', 'GSWP3_W5E5']
pcps_all = ['pcp2.5', 'pcp1.6','2.6', 'via_winterprcp']
for qc in ['qc0', 'qc3']:
for pcp, clim in zip(pcps_all, clims_all):
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
if qc == 'qc0':
border = 80
alpha_min = 0
if 'W5E5' in clim:
border = 160
else:
border = 160
alpha_min = 0.7
if '2.6' == pcp:
climi = clim +'_cte'
elif 'via_winterprcp' == pcp:
climi = clim + '_var'
else:
climi = clim
try:
if 'match_geod_pergla' in match and qc=='qc3':
pass
else:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdt_OGGM_{dict_key_short}', marker=markers[climi],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp]-alpha_min,
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
dfh.plot(ax=ax, y='dmdt', yerr='err_dmdt', marker='s', linestyle='none', markersize=12,
label = 'geodetic_observation_Hugonnet2021', color = 'orange', alpha = 0.5)
dfh.plot(ax=ax, y='dmdt_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation (only common running glaciers)', color = 'violet', alpha = 1, markeredgewidth=5)
plt.ylabel(r'regional total mass change (dmdt, Gt year$^{-1}$)')
plt.xlabel('RGI region')
plt.xticks(np.arange(0,18,1), dfh.index.values)
f = ax.get_legend()
f.remove()
###
ax = plt.subplot2grid((n, n), (9, 0), colspan=8, rowspan=8)
#f, ax = plt.subplots()
dfh.plot(ax=ax, y='dmdtda', yerr='err_dmdtda', marker='s', linestyle='none', markersize=12,
label = 'geodetic observation (Hugonnet et al., 2021)', color = 'orange', alpha = 0.7)
dfh.plot(ax=ax, y='dmdtda_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation (only common running glaciers)', color = 'violet', alpha = 1, markeredgewidth=5)
for qc in ['qc0', 'qc3']:
for pcp, clim in zip(pcps_all, clims_all):
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
if qc == 'qc0':
border = 80
alpha_min = 0
if 'W5E5' in clim:
border = 160
else:
border = 160
alpha_min = 0.7
if '2.6' == pcp:
climi = clim +'_cte'
elif 'via_winterprcp' == pcp:
climi = clim + '_var'
else:
climi = clim
try:
if 'match_geod_pergla' in match and qc=='qc3':
pass
else:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdtda_OGGM_{dict_key_short}', marker=markers[climi],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp]-alpha_min,
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
#t=ax.get_legend_handles_labels()
#t[0] = [t[0][-1]].append(t[:-1])
plt.ylabel(r'area-weighted specific mass balance (dmdtda, kg m$^{-2}$ year$^{-1}$)')
plt.xlabel('RGI region')
plt.xticks(np.arange(0, 18, 1), dfh.index.values)
f = ax.get_legend()
f.remove()
plt.title('area-weighted mean MB from "fixed_geometry_mass_balance"')
##
#plt.figure(figsize=(20,10))
ax = plt.subplot2grid((n, n), (9, 9), colspan=8, rowspan=8)
for qc in ['qc0', 'qc3']:
for pcp, clim in zip(pcps_all, clims_all):
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
if qc == 'qc0':
border = 80
alpha_min = 0
if 'W5E5' in clim:
border = 160
else:
border = 160
alpha_min = 0.7
if '2.6' == pcp:
climi = clim +'_cte'
elif 'via_winterprcp' == pcp:
climi = clim + '_var'
else:
climi = clim
try:
if 'match_geod_pergla' in match and qc=='qc3':
pass
else:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdtda_dyna_OGGM_{dict_key_short}', marker=markers[climi],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp]-alpha_min,
label=dict_key_short)
except:
pass
dfh.plot(ax=ax, y='dmdtda', yerr='err_dmdtda', marker='s', linestyle='none', markersize=12,
label = 'geodetic observation (Hugonnet et al., 2021)', color = 'orange', alpha = 0.7)
dfh.plot(ax=ax, y='dmdtda_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation (only common running glaciers)', color = 'violet', alpha = 1, markeredgewidth=5)
plt.ylabel(r'dynamic specific mass balance (dmdtda, kg m$^{-2}$ year$^{-1}$)')
plt.xlabel('RGI region')
plt.title('using volume changes from "historical_run_output_extended"')
plt.xticks(np.arange(0,18,1), dfh.index.values);
plt.ylim([-1600,-50])
ax.legend(framealpha=0.5, ncol=1, loc=3, bbox_to_anchor=(0.17,1.3),
title=f'mean of 2000-2020: qc0 vs qc3 (only glaciers that are\nrunning for all preprocessed gdirs used, i.e. {len(all_running_rgis)*100/len(pd_working):0.1f}%)')
#######
ax = plt.subplot2grid((n, n), (1, 9), rowspan=6, colspan=1)
rgi_index = ['{:02d}'.format(int(rgi_reg)) for rgi_reg in np.arange(1,19)]
dfh.loc['all_without_19', 'dmdt'] = dfh.loc[rgi_index].dmdt.sum()
dfh.loc['all_without_19', 'err_dmdt'] = dfh.loc[rgi_index].err_dmdt.sum()
dfh.loc['all_without_19', 'period'] = dfh.period[0]
for qc in ['qc0', 'qc3']:
for pcp, clim in zip(pcps_all, clims_all):
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
if qc == 'qc0':
border = 80
if 'W5E5' in clim:
border = 160
else:
border = 160
if '2.6' == pcp:
climi = clim +'_cte'
elif 'via_winterprcp' == pcp:
climi = clim + '_var'
else:
climi = clim
try:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}'
dfh.loc['all_without_19', f'dmdt_OGGM_{dict_key_short}'] = dfh.loc[rgi_index, f'dmdt_OGGM_{dict_key_short}'].sum()
except:
pass
###
ax = plt.gca()
#ax.errorbar(['geodetic observation (Hugonnet et al., 2021)'],
# df_dmdt_all.dmdt, yerr=df_dmdt_all.err_dmdt,
# marker='s', color='orange', alpha=0.8)
df_dmdt_all = dfh.loc['all_without_19'][1:].dropna()
plt.axhline(df_dmdt_all.dmdt, color='orange', alpha = 0.3)
plt.axhspan(df_dmdt_all.dmdt-df_dmdt_all.err_dmdt,
df_dmdt_all.dmdt + df_dmdt_all.err_dmdt, alpha = 0.1, color='orange')
plt.axhline((pd_geodetic_running.dmdtda*1e3*pd_geodetic_running.area).sum()*1e-12, color='violet', alpha=0.5, linewidth=5)
for qc in ['qc0', 'qc3']:
for pcp, clim in zip(pcps_all, clims_all):
for match in ['no_match', 'match_geod', 'match_geod_pergla', 'match_geod_pergla_massredis']: #
if qc == 'qc0':
border = 80
alpha_min = 0
if 'W5E5' in clim:
border = 160
else:
border = 160
alpha_min = 0.7
if '2.6' == pcp:
climi = clim +'_cte'
elif 'via_winterprcp' == pcp:
climi = clim + '_var'
else:
climi = clim
try:
if 'match_geod_pergla' in match and qc=='qc3':
pass
else:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
ax.errorbar(x=[dict_key_short], # [0]
y=df_dmdt_all[f'dmdt_OGGM_{dict_key_short}'],
marker=markers[climi],
linestyle='none', markersize=8, color=colors[match], alpha=alphas[exp]-alpha_min,
label=dict_key_short)
except:
pass
plt.ylabel(r'world-wide total mass-change (dmdt, Gt year$^{-1}$)')
plt.xticks(ticks=[3.5],labels=['all without RGI 19'])
plt.xlim([-1,9.5])
plt.grid(axis='x')
plt.savefig(f'dmdtda_dmdt_for_prepro_level_gdirs_qc0_vs_qc3.png', bbox_inches='tight',
pad_inches=0)
plt.savefig(f'dmdtda_dmdt_for_prepro_level_gdirs_qc0_vs_qc3.pdf', bbox_inches='tight',
pad_inches=0.2)
plt.figure(figsize=(10,18))
ax = plt.gca()
#f, ax = plt.subplots()
dfh.plot(ax=ax, y='dmdtda', yerr='err_dmdtda', marker='s', linestyle='none', markersize=12,
label = 'geodetic observation (Hugonnet et al., 2021)', color = 'orange', alpha = 0.7)
dfh.plot(ax=ax, y='dmdtda_geodetic_only_running_glaciers', marker="_", linestyle='none', markersize=30,
label = 'geodetic observation (only common running glaciers)', color = 'violet', alpha = 1, markeredgewidth=5)
for qc in ['qc0']: #, 'qc3']:
for pcp, clim in zip(pcps_all, clims_all):
for match in ['match_geod_pergla', 'match_geod_pergla_massredis']: #
if qc == 'qc0':
border = 80
alpha_min = 0
if 'W5E5' in clim:
border = 160
else:
border = 160
alpha_min = 0.7
if '2.6' == pcp:
climi = clim +'_cte'
elif 'via_winterprcp' == pcp:
climi = clim + '_var'
else:
climi = clim
try:
if 'match_geod_pergla' in match and qc=='qc3':
pass
else:
dict_key_short = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}' #_rgi_{rgi_reg}'
dfh.plot(ax=ax, y=f'dmdtda_OGGM_{dict_key_short}', marker=markers[climi],
linestyle='none', markersize=8, color = colors[match], alpha = alphas[exp]-alpha_min,
label=f'{exp}_{clim}_{pcp}_{match}')
except:
pass
#t=ax.get_legend_handles_labels()
#t[0] = [t[0][-1]].append(t[:-1])
plt.ylabel(r'area-weighted specific mass balance (dmdtda, kg m$^{-2}$ year$^{-1}$)')
plt.xlabel('RGI region')
plt.xticks(np.arange(0, 18, 1), dfh.index.values[:-1])
plt.legend()
#f.remove()
plt.title('area-weighted mean MB from "fixed_geometry_mass_balance"')
/home/users/lschuster/.local/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1061: UserWarning: Warning: converting a masked element to nan. x = np.asanyarray(x)
Text(0.5, 1.0, 'area-weighted mean MB from "fixed_geometry_mass_balance"')
rgi_reg='11'
rgi_reg_int = 11
plt.figure(figsize=(20,20))
#for rgi_reg_int, rgi_reg in enumerate(np.arange(1,):
#plt.subplot(6,2,rgi_reg_int)
plt.subplot(221)
fd = '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod_pergla/RGI62/b_080/L5/summary/'
pd_test = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == rgi_reg_int]
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
plt.boxplot(pd_test.loc[2000:,all_running_rgis_reg].mean()-pd_geodetic_running_reg.loc[all_running_rgis_reg].dmdtda *1e3)
#plt.xlabel('oggm SMB match_geod')
#plt.ylabel('geodetic SMB')
plt.ylabel('oggm match_geod - geodetic SMB')
plt.title('old: {}'.format(rgi_reg))
plt.subplot(222)
fd = '/home/www/oggm/gdirs/oggm_v1.6/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod_pergla/RGI62/b_080/L5/summary/'
pd_test = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == rgi_reg_int]
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
plt.boxplot(pd_test.loc[2000:,all_running_rgis_reg].mean()-pd_geodetic_running_reg.loc[all_running_rgis_reg].dmdtda *1e3)
#plt.xlabel('oggm SMB match_geod')
#plt.ylabel('geodetic SMB')
plt.ylabel('oggm match_geod - geodetic SMB')
plt.title('new: {}'.format(rgi_reg))
#
plt.subplot(223)
rgi_reg='05'
rgi_reg_int = 5
#for rgi_reg_int, rgi_reg in enumerate(np.arange(1,):
#plt.subplot(6,2,rgi_reg_int)
fd = '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod_pergla/RGI62/b_080/L5/summary/'
pd_test = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == rgi_reg_int]
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
plt.boxplot(pd_test.loc[2000:,all_running_rgis_reg].mean()-pd_geodetic_running_reg.loc[all_running_rgis_reg].dmdtda *1e3)
#plt.xlabel('oggm SMB match_geod')
#plt.ylabel('geodetic SMB')
plt.ylabel('oggm match_geod - geodetic SMB')
plt.title('old: {}'.format(rgi_reg))
#
plt.subplot(224)
rgi_reg='05'
rgi_reg_int = 5
#for rgi_reg_int, rgi_reg in enumerate(np.arange(1,):
#plt.subplot(6,2,rgi_reg_int)
fd = '/home/www/oggm/gdirs/oggm_v1.6/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod_pergla/RGI62/b_080/L5/summary/'
pd_test = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == rgi_reg_int]
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
plt.boxplot(pd_test.loc[2000:,all_running_rgis_reg].mean()-pd_geodetic_running_reg.loc[all_running_rgis_reg].dmdtda *1e3)
#plt.xlabel('oggm SMB match_geod')
#plt.ylabel('geodetic SMB')
plt.ylabel('oggm match_geod - geodetic SMB')
plt.title('new: {}'.format(rgi_reg))
Text(0.5, 1.0, 'new: 05')
plt.figure(figsize=(20,10))
plt.subplot(121)
rgi_reg='05'
rgi_reg_int = 5
#for rgi_reg_int, rgi_reg in enumerate(np.arange(1,):
#plt.subplot(6,2,rgi_reg_int)
fd = '/home/www/fmaussion/runs/new_gdirs/oggm_v1.6/L3-L5_files/GSWP3_W5E5/elev_bands/qc0/match_geod_pergla_2.6/RGI62/b_160/L5/summary/'
pd_test = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == rgi_reg_int]
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
plt.boxplot(pd_test.loc[2000:,all_running_rgis_reg].mean()-pd_geodetic_running_reg.loc[all_running_rgis_reg].dmdtda *1e3)
#plt.xlabel('oggm SMB match_geod')
#plt.ylabel('geodetic SMB')
plt.ylabel('oggm match_geod - geodetic SMB')
plt.title('W5E5: cte prcp. fac {}'.format(rgi_reg))
plt.subplot(122)
rgi_reg='05'
rgi_reg_int = 5
#for rgi_reg_int, rgi_reg in enumerate(np.arange(1,):
#plt.subplot(6,2,rgi_reg_int)
fd = '/home/www/fmaussion/runs/new_gdirs/oggm_v1.6/L3-L5_files/GSWP3_W5E5/elev_bands/qc0/match_geod_pergla_winterprcp/RGI62/b_160/L5/summary/'
pd_test = pd.read_csv(fd + f'fixed_geometry_mass_balance_{rgi_reg}.csv'.format(rgi_reg), index_col=0, low_memory=False)
pd_geodetic_running_reg = pd_geodetic_running[pd_geodetic_running.reg == rgi_reg_int]
all_running_rgis_reg = pd_working.loc[pd_working.rgi_reg==rgi_reg_int]['all_running_rgis'].dropna().index
plt.boxplot(pd_test.loc[2000:,all_running_rgis_reg].mean()-pd_geodetic_running_reg.loc[all_running_rgis_reg].dmdtda *1e3)
#plt.xlabel('oggm SMB match_geod')
#plt.ylabel('geodetic SMB')
plt.ylabel('oggm match_geod - geodetic SMB')
plt.title('W5E5: variable prcp. fac {}'.format(rgi_reg))
Text(0.5, 1.0, 'W5E5: variable prcp. fac 05')
match_geod_pergla
calibration!¶rgi_reg = 'all_without_19'
pcps = ['pcp2.5']
clims = ['CRU']
dvol_wrong = {}
print('missing combination options:')
for qc in ['qc0']:
for exp in ['elev_bands']: #, 'centerlines']:
for match in ['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.4/'
#else:
# take the wrong preprocessed gdirs just a s a comparison
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_wrong[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!!!)
except:
print(dict_key)
missing combination options:
plt.rc('font', size=18)
plt.figure(figsize=(26,14))
alpha = 1
for match in ['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])
label = f'{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_wrong'
plt.plot(dvol_wrong[dict_key].time[dvol[dict_key].time>=2000], dvol_wrong[dict_key][dvol[dict_key].time>=2000].values/1e9,
label=label,alpha=alpha-0.5,
ls=ls, color=colors[match])
except:
print('missing: {}'.format(dict_key))
plt.legend(title = 'new vs wrong calibration')
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);