#!/usr/bin/env python
# coding: utf-8
# # "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](https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~lschuster/error_analysis/error_analysis_v1.ipynb?flush_cache=true#Analysis-for-Level-5-pre-processing-directories!).
#
# 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:
# - [historical volume changes (2000-2020) for qc3 (without match_geod_pergla)](#id-volume-evolution-qc3)
# - [comparison of historical volume changes (2000-2020): qc0 vs qc3](#id-volume-evolution-qc3-vs-qc0)
# - [comparison of historial volume changes: looking at all glaciers vs only the "common running" glaciers](#id-volume-evolution-all-glaciers-vs-common-running-glaciers)
# - [regional mass change or specific mass change with qc3](#id-dmdt-dmdtda-qc3)
# - [regional mass change or specific mass change: comparison qc0 vs qc3](#id-dmdt-dmdtda-qc3-vs-qc0)
#
# 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.
# - [differences in historical volume projections between the new and the old wrong preprocessed gdirs from match_geod_perglac](#New-vs-old-wrong-match_geod_pergla-calibration!)
# 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](#id-dmdt-dmdtda-qc3-vs-qc0))
# - in the relative volume change: for no_match (qc0 has a much more negative volume change than qc3, see [dmdt plot below](#id-dmdt-dmdtda-qc3-vs-qc0)!)
#
# 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)](https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~lschuster/error_analysis/vas_error_analysis.ipynb?flush_cache=true#id-volume-proj-single-gcm)!
# 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);
#
# - as there are more glaciers that run without errors with ERA5 than with CRU, we have in total a larger volume for ERA5 compared to CRU
# - only when using match_geod_pergla, the difference between ERA5 and CRU are small. This is because the error amount difference between CRU and ERA5 for this match option is very small!
# - see more in [this discussion and plot of the error analysis notebook](https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~lschuster/error_analysis/error_analysis_v1.ipynb?flush_cache=true#id-total-error-area-summary)
# In[ ]:
# ### 1.4 Compare climate datasets (CRU, ERA5, W5E5 with constant & variable prcp. fac)
# In[12]:
# 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)
# In[13]:
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)
# ## 2. Regional total mass change and specific mass balance using "fixed_geometry_mass_balance" and average "historical_run_output_extended"
# In[60]:
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')
# In[61]:
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]
# In[62]:
dfh.columns
# In[42]:
# In[63]:
# 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
# In[64]:
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)
# ### 2.1 only look at dmdtda and at dmdt with qc3 (climate quality check and correction if necessary)
# - only for no_match and match_geod_pergla, without W5E5
# In[69]:
dfh['dmdtda_OGGM_elev_bands_pcp1.6_ERA5_match_geod_pergla_qc0_b80']
# In[72]:
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
# - when looking at all glaciers that run for one preprocessed gdirs option and using match_geod, the non-dynamical SMB estimates of OGGM will match perfectly to the regional geodetic calibration data of Hugonnet (see: [the fixed-area SMB section of the "compare_with_geod" notebook](https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/compare_with_geod.ipynb?flush_cache=true))
#
# - the dynamic SMB matches less good (specifically visible for match_geod) because the changing area was neglected during the MB calibration.
# - the single influence of that can be seen by comparing always only one gdirs option with all running glaciers to match_geod (as above, see [the dynamic SMB section of the "compare_with_geod" notebook](https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/compare_with_geod.ipynb?flush_cache=true#Now-compare-the-SMB-of-the-dynamics-run,-which-is-not-perfect-anymore))
#
# Somehow naturally, the **no_match** option is most different to the geodetic estimates (as here only the direct glaciological reference data are used):
# - no_match climate differences neglecting the dynamic influence:
# - RGI regions 01, 08, 13-17: ERA5 has a less negative specific MB than CRU, for others CRU
# - RGI regions 04, 05, 07, 09-12, 18: ERA5 has a more negative specific MB than CRU, for others CRU
# - no_match climate differences including the dynamic influence:
# - similar as when neglecting dynamic influence, however differences are less strong!
#
# 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**
# - in terms of mass change, of course, using the geodetic estimates for only the common running glaciers gives a smaller total mass loss
# - for dmdtda, for most RGI regions, using **geodetic estimates for only the common running glaciers** results in a rather slightly less negative dmdtda
# ### 2.2 compare dmdtda and dmdt between qc0 (no climate check and correction) and qc3 (climate quality check and correction if necessary) for elev_bands only!
# In[ ]:
# In[ ]:
# In[76]:
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)
#
# - no_match with qc0 overestimates the total world-wide mass change in the period 2000 to 2020 even more than no_match with qc3!!!
# - in all other options, qc0 is similar to qc3, but always result globally in a slightly smaller mass loss!!!
# - when using **match_geod_pergla**, every glacier that runs should perfectly match the dmdtda of that glacier for the non-dynamic SMB!
# - we will check this further in the following plots:
# In[77]:
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"')
# - this is the same plot as before, but now only with the match_geod_pergla gdirs that we want to check, they seem to match perfectly over the mean!
# - below is also a check if they match on the per-glacier level!
# In[78]:
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))
# - the problem from the old preprocessed gdirs from oggm_v1.4 is solved in the new preprocessed gdirs (oggm_v1.6)!
# In[ ]:
# In[85]:
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))
# In[ ]:
# # Appendix
# ## New vs old wrong `match_geod_pergla` calibration!
# - How different is it for the historical projections?is
#
# In[25]:
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)
# In[26]:
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);
# In[ ]: