Convert "raw" Hugonnet et al (2021) geodetic per glacier dataset to a corrected dataset without outliers and missing glaciers¶

df

len(dfs)

not_rgi = np.array(['RGI60' not in s for s in dfs.index])
len(dfs.loc[not_rgi]), len(dfs.loc[~not_rgi])

not_rgi = np.array(['RGI60' not in s for s in df_new.index])
df_new.loc[not_rgi]

df_new.loc[not_rgi].reg.unique()

template

n = len(df_new.loc[rgi['O1Region']!='12'][df_new.loc[rgi['O1Region']!='12'].period == '2000-01-01_2020-01-01'])
print('glaciers outside RGI region 12 inside geodetic dataset : {}'.format(n))
print('glaciers outside RGI region 12 in total from rgi : {}'.format(len(rgi.loc[rgi['O1Region']!='12'])))
print('glaciers in total in rgi : {}'.format(len(rgi)))

df_new = df_new.rename_axis('rgiid').sort_values(by = ['rgiid', 'period'])
df_new

# there are some glaciers on connectivity level 2
len(rgi.loc[rgi['Connect'] == 2])

# check: they are only in RGI region 05 (Greenland)
rgi.loc[rgi['Connect'] == 2]['O1Region'].unique()

len(rgi)

dfs.err_dmdtda.quantile([0.001, 0.01, 0.5, 0.99, 0.999])

dfs.err_dmdtda.max()

# SHOULD IT BE AREA WEIGHTED???? Maybe not
all_sigma_mean = dfs['err_dmdtda'].mean()
all_sigma_std = dfs['err_dmdtda'].std()
all_sigma_threshold = all_sigma_mean + 3 * all_sigma_std
print('all sigma threshold (global estimate):', np.round(all_sigma_threshold, 2))

periods = ['2000-01-01_2010-01-01', '2000-01-01_2020-01-01', '2010-01-01_2020-01-01']
# We filter and fill all periods based on the 20 year threshold
df_filled = df_new.copy()[['period', 'area', 'dmdtda', 'err_dmdtda', 'reg']]
df_filled['is_cor'] = False
for reg in range(1, 20):
    
    dfs_subset = dfs.loc[dfs['reg'] == reg]
    
    # Too high of sigma causes large issues for model
    # compute regional threshold for every region
    reg_sigma_mean = dfs_subset['err_dmdtda'].mean()
    reg_sigma_std = dfs_subset['err_dmdtda'].std()
    reg_sigma_threshold = reg_sigma_mean + 3 * reg_sigma_std
    
    # Don’t penalize regions that are well-measured, so use all threshold as minimum:
    # if the regional threshold is smaller than the global threshold,
    # we use the global threshold, otherwise we use the regional threshold
    if reg_sigma_threshold < all_sigma_threshold:
        reg_sigma_threshold = all_sigma_threshold
        
    to_replace = dfs_subset.loc[dfs_subset['err_dmdtda'] > reg_sigma_threshold]
    to_keep = dfs_subset.loc[dfs_subset['err_dmdtda'] <= reg_sigma_threshold]
    
    print(reg, len(dfs_subset), len(to_replace),  np.round(reg_sigma_threshold, 2))
        
    df_filled.loc[to_replace.index, 'dmdtda'] = np.NaN
    df_filled.loc[to_replace.index, 'err_dmdtda'] = np.NaN
    
    # Replace nan values - SHOULD BE AREA WEIGHTED????
    for period in periods:
        
        # indices of glaciers without dmdtda data (because missing or outlier)
        loc_no = (df_filled['reg'] == reg) & (df_filled['dmdtda'].isnull()) & (df_filled['period'] == period)
        # indices of glaciers with dmdtda data
        loc_yes = (df_filled['reg'] == reg) & (~ df_filled['dmdtda'].isnull()) & (df_filled['period'] == period)
        
        # replace the nan-values from the missing and outlier glaciers with the regional mean and standard deviation mean 
        # Note that every glacier without "usable" geodetic data will get the same dmdtda and err_dmdtda 
        # (independent of the glacier characteristics such as size, elevation) !!!
        df_filled.loc[loc_no,  'dmdtda'] = df_filled.loc[loc_yes,  'dmdtda'].mean()
        df_filled.loc[loc_no,  'err_dmdtda'] = df_filled.loc[loc_yes,  'err_dmdtda'].mean()
        df_filled.loc[loc_no,  'is_cor'] = True

dfs_filled.err_dmdtda.quantile([0.001, 0.01, 0.5, 0.99, 0.999])

dfs_filled.err_dmdtda.max()

dfs_filled.dmdtda.max()

fig, ax = plt.subplots()
dfs_filled.dmdtda.plot(kind='hist', bins=100, bottom=0.1);
ax.set_yscale('log')
plt.title(f'Distribution of the specific mass balance after correction for all glaciers (n={len(dfs_filled)})')
plt.xlabel(r'specific mass balance (kg m$^{-2}$)');

# relative amount of glacier area with a positive specific mass balance
dfs_filled.loc[dfs_filled.dmdtda > 0].area.sum() / dfs_filled.area.sum()

fig, ax = plt.subplots(ncols=3, figsize=(18,6))
plt.suptitle(f'Distribution of the specific mass balance after correction for all glaciers (n={len(dfs_filled)})')
for k,period in enumerate([periods[1], periods[0], periods[2]]):
    for _,period2 in enumerate([periods[1], periods[0], periods[2]]):
        dfs_filled = df_filled.loc[df_new.period == period2]
        dfs_filled.dmdtda.plot(kind='hist', bins=100, bottom=0.1, alpha=0.3, ax = ax[k], color = 'grey', label='other periods');
    dfs_filled = df_filled.loc[df_new.period == period]
    dfs_filled.dmdtda.plot(kind='hist', bins=100, bottom=0.1, alpha=1, ax = ax[k], label = 'actual period');
    ax[k].set_yscale('log')
    ax[k].set_title('period: {},\n max error: {:0.2f}, max dmdtda: {:0.2f},\n glacier area with dmdtda above zero : {:0.2f}%'.format(period,
                                                                                                          dfs_filled.err_dmdtda.max(),
                                                                                                          dfs_filled.dmdtda.max(),
                                                                                                          dfs_filled.loc[dfs_filled.dmdtda > 0].area.sum()*100 / dfs_filled.area.sum()))
    ax[k].legend()
    ax[k].set_xlabel(r'specific mass balance (kg m$^{-2}$)');
plt.tight_layout()

# check: are there any glaciers left with nan-values?
df_filled.loc[df_filled['dmdtda'].isnull()]
# No!

df_filled_per = df_filled[df_filled['period']=='2000-01-01_2020-01-01']
n_cor = len(df_filled_per.loc[df_filled_per['is_cor']])
n_total = len(df_filled_per)
area_cor = df_filled_per.loc[df_filled_per['is_cor']].area
area_cor_sum = area_cor.sum()
area_cor_med = area_cor.median()
print(f'{n_cor} out of {n_total} glaciers were filled up with regional mean data, hence {np.round(n_cor*100/n_total,2)}%.')
print(f'{np.round(area_cor_sum*100/df_filled_per.area.sum(),2)}% of the glacier area was refilled with regional mean data')

# same but using instead regional medians for corrections
df_filled_per_med = df_filled_med[df_filled_med['period']=='2000-01-01_2020-01-01']

plt.figure(figsize=(20,10))
plt.subplot(121)
plt.plot(df_filled_per.loc[~df_filled_per['is_cor']].area,df_filled_per.loc[~df_filled_per['is_cor']].dmdtda, '.', label='no correction necessary')
plt.plot(df_filled_per.loc[df_filled_per['is_cor']].area,df_filled_per.loc[df_filled_per['is_cor']].dmdtda, '.', 
         label='corrected by regional means')
plt.plot(df_filled_per_med.loc[df_filled_per_med['is_cor']].area,
         df_filled_per_med.loc[df_filled_per['is_cor']].dmdtda, 'x',alpha=0.4, ms=4,
         label='corrected by regional medians')
plt.legend()
plt.xlabel('area (m2)')
plt.ylabel('dmdtda')

plt.subplot(122)
plt.plot(df_filled_per.loc[~df_filled_per['is_cor']].area,df_filled_per.loc[~df_filled_per['is_cor']].err_dmdtda, '.', label='no correction necessary')
plt.plot(df_filled_per.loc[df_filled_per['is_cor']].area,df_filled_per.loc[df_filled_per['is_cor']].err_dmdtda, '.', 
         label='corrected by regional means')
plt.plot(df_filled_per_med.loc[df_filled_per_med['is_cor']].area,df_filled_per_med.loc[df_filled_per_med['is_cor']].err_dmdtda,
         'x', alpha=0.4, ms=4,label='corrected by regional medians')
plt.legend()
plt.xlabel('area (m2)')
plt.ylabel('err_dmdtda');

plt.figure(figsize=(15,3))
plt.plot(df_new_per.area, df_new_per.dhdt, '.')
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area,0.025), ls='--')
plt.axvline(df_filled[df_filled.is_cor].area.median())
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area,0.975), ls='--')
plt.xlabel('area (m2)')
plt.ylabel('dhdt')

plt.figure(figsize=(15,10))
plt.plot(df_new_per_small.area, df_new_per_small.dhdt, '.')
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area,0.025), ls='--')
plt.axvline(df_filled[df_filled.is_cor].area.median())
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area,0.975), ls='--')
plt.xlabel('area (m2)')
plt.ylabel('dhdt')

plt.figure(figsize=(15,10))
plt.plot(df_new_per_small.area, df_new_per_small.err_dhdt, '.')
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area,0.025), ls='--')
plt.axvline(df_filled[df_filled.is_cor].area.median())
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area,0.975), ls='--')
plt.xlabel('area (m2)')
plt.ylabel('dhdt')

plt.figure(figsize=(15,3))
plt.plot(df_new_per.area, df_new_per.dvoldt, '.')
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area, 0.025), ls='--')
plt.axvline(df_filled[df_filled.is_cor].area.median())
plt.axvline(np.quantile(df_filled[~df_filled.is_cor].area, 0.975), ls='--')
plt.xlabel('area (m2)')
plt.ylabel('dvoldt (m3)')

plt.figure(figsize=(15,10))
plt.subplot(121)
reg = LinearRegression().fit(np.atleast_2d(df_new_per_small.area.values).T, df_new_per_small.dvoldt)
missing_area_glaciers = np.atleast_2d(df_filled_per[df_filled_per.is_cor].area.values).T
dvoldt_pred_cor = reg.predict(missing_area_glaciers)
dvoldt_pred_test = reg.predict(np.atleast_2d(df_new_per.loc[df_new_per_small.index].area.values).T)
plt.plot(df_new_per_small.area, df_new_per_small.dvoldt, '.', label='world-wide RGI glaciers')
plt.plot(missing_area_glaciers, dvoldt_pred_cor, 'o', label='filled up missing glaciers')
plt.axvline(np.quantile(df_filled_per[df_filled_per.is_cor].area,0.0025), ls='--', label='00.25%-area of missing glaciers')
plt.axvline(df_filled_per[df_filled_per.is_cor].area.median(), label='median of missing glaciers')
plt.axvline(np.quantile(df_filled_per[df_filled_per.is_cor].area,0.9975), ls='--', label='99.75%-area of missing glaciers')
plt.xlabel('area (m2)')
plt.ylabel('dvoldt')
plt.plot(df_new_per.loc[df_new_per_small.index].area.values, dvoldt_pred_test, '.', color = 'grey', ms=3, alpha = 0.2, label='fitted')
plt.title('R2-value: {:0.2f}'.format(sklearn.metrics.r2_score(df_new_per_small.dvoldt,
                               dvoldt_pred_test)))
plt.legend()

plt.subplot(122)
reg2 = LinearRegression().fit(np.atleast_2d(df_new_per_small.area.values).T, df_new_per_small.err_dvoldt)
missing_area_glaciers = np.atleast_2d(df_filled_per[df_filled_per.is_cor].area.values).T
err_dvoldt_pred_cor = reg2.predict(missing_area_glaciers)
err_dvoldt_pred_test = reg2.predict(np.atleast_2d(df_new_per.loc[df_new_per_small.index].area.values).T)


plt.plot(df_new_per_small.area, df_new_per_small.err_dvoldt, '.', label='world-wide RGI glaciers')
plt.plot(missing_area_glaciers, err_dvoldt_pred_cor, 'o', label='filled up missing glaciers')
plt.axvline(np.quantile(df_filled_per[df_filled_per.is_cor].area,0.0025), ls='--', label='00.25%-area of missing glaciers')
plt.axvline(df_filled_per[df_filled_per.is_cor].area.median(), label='median of missing glaciers')
plt.axvline(np.quantile(df_filled_per[df_filled_per.is_cor].area,0.9975), ls='--', label='99.75%-area of missing glaciers')
plt.xlabel('area (m2)')
plt.ylabel('err_dvoldt')
plt.plot(df_new_per.loc[df_new_per_small.index].area.values, err_dvoldt_pred_test, '.', color = 'grey', ms=3, alpha = 0.2, label='fitted')
plt.title('R2-value: {:0.2f}'.format(sklearn.metrics.r2_score(df_new_per_small.err_dvoldt,
                               err_dvoldt_pred_test)))
plt.legend()
#plt.savefig('dvoldt_missing_fit_global.png')