import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from oggm import utils
import xarray as xr
import glob, os
dfz = pd.read_csv(utils.get_demo_file('zemp_ref_2006_2016.csv'), index_col=0)
dfh = pd.read_csv(utils.get_demo_file('table_hugonnet_regions_10yr_20yr_ar6period.csv'))
dfh = dfh.loc[dfh.period == '2000-01-01_2020-01-01'].set_index('reg')
dfz['SMB_ZEMP'] = dfz['SMB'] * 1000
dfz['SMB_ZEMP_err'] = dfz['SMB_err'] * 1000
dfz['SMB_HUGUONNET'] = dfh['dmdtda']
dfz['SMB_HUGUONNET_err'] = dfh['err_dmdtda']
dfz.index = ['{:02d}'.format(rgi_reg) for rgi_reg in dfz.index]
f, ax = plt.subplots()
dfz.plot(ax=ax, y='SMB_ZEMP', yerr='SMB_ZEMP_err', marker='o', linestyle='none');
dfz.plot(ax=ax, y='SMB_HUGUONNET', yerr='SMB_HUGUONNET_err', marker='o', linestyle='none');
idir = '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/'
exps = glob.glob(idir + '*/*/qc*/pcp*/match_geod')
exps
['/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/ERA5/elev_bands/qc3/pcp1.6/match_geod', '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/ERA5/elev_bands/qc3/pcp1.8/match_geod', '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc3/pcp2.5/match_geod', '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod', '/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CERA+ERA5/elev_bands/qc3/pcp1.6/match_geod']
for exp in exps:
dfz['SMB_OGGM'] = np.NaN
dfz['SMB_OGGM_calv'] = np.NaN
odir = os.path.join(exp, 'RGI62', 'b_080', 'L5', 'summary')
title = exp.replace(idir, '').replace('/match_geod', '').replace('/elev_bands/', '/').replace('/', ' ')
for rgi_reg in np.arange(1, 20):
rgi_reg = '{:02d}'.format(rgi_reg)
try:
df = pd.read_csv(odir + '/fixed_geometry_mass_balance_{}.csv'.format(rgi_reg), index_col=0, low_memory=False)
dfs = pd.read_csv(odir + '/glacier_statistics_{}.csv'.format(rgi_reg), index_col=0, low_memory=False)
except FileNotFoundError:
print(odir, rgi_reg, 'MISSING')
continue
df = df.dropna(axis=0, how='all')
df = df.dropna(axis=1, how='all')
odf = pd.DataFrame(df.loc[2000:].mean(), columns=['SMB'])
odf['AREA'] = dfs.rgi_area_km2
dfz.loc[rgi_reg, 'SMB_OGGM'] = np.average(odf['SMB'], weights=odf['AREA'])
f, ax = plt.subplots()
dfz.plot(ax=ax, y='SMB_OGGM', marker='o', linestyle='none', markersize=9);
dfz.plot(ax=ax, y='SMB_HUGUONNET', yerr='SMB_HUGUONNET_err', marker='o', linestyle='none');
plt.title(title);
/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc3/pcp2.5/match_geod/RGI62/b_080/L5/summary 19 MISSING /home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod/RGI62/b_080/L5/summary 19 MISSING
for exp in exps:
odir = os.path.join(exp, 'RGI62', 'b_080', 'L5', 'summary')
title = exp.replace(idir, '').replace('/match_geod', '').replace('/elev_bands/', '/').replace('/', ' ')
odfs = pd.DataFrame()
for rgi_reg in np.arange(1, 20):
rgi_reg = '{:02d}'.format(rgi_reg)
try:
dfs = pd.read_csv(odir + '/glacier_statistics_{}.csv'.format(rgi_reg), index_col=0, low_memory=False)
except FileNotFoundError:
print(odir, rgi_reg, 'MISSING')
continue
dfs = dfs[['mb_bias_before_geodetic_corr', 'mb_bias', 'rgi_area_km2']].dropna()
odfs.loc[rgi_reg, 'avg_bias_before_corr'] = np.average(dfs['mb_bias_before_geodetic_corr'], weights=dfs['rgi_area_km2'])
odfs.loc[rgi_reg, 'avg_bias_after_corr'] = np.average(dfs['mb_bias'], weights=dfs['rgi_area_km2'])
odfs.loc[rgi_reg, 'bias_corr'] = (dfs.mb_bias.iloc[0] - dfs.mb_bias_before_geodetic_corr.iloc[0])
odfs.loc[rgi_reg, 'area'] = dfs.rgi_area_km2.sum()
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
odfs.plot(ax=ax1, y='bias_corr', marker='o', linestyle='none', markersize=9);
ax1.set_title(title + ' (avg corr: {:.0f})'.format(np.average(odfs.bias_corr, weights=odfs['area'])));
odfs.plot(ax=ax2, y='avg_bias_before_corr', marker='o', linestyle='none', markersize=9);
odfs.plot(ax=ax2, y='avg_bias_after_corr', marker='o', linestyle='none', markersize=9);
ax2.set_title('avg abs bias before: {:.0f} after: {:.0f}'.format(np.average(odfs.avg_bias_before_corr.abs(), weights=odfs.area),
np.average(odfs.avg_bias_after_corr.abs(), weights=odfs.area)));
/home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc3/pcp2.5/match_geod/RGI62/b_080/L5/summary 19 MISSING /home/www/oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc0/pcp2.5/match_geod/RGI62/b_080/L5/summary 19 MISSING
for exp in exps:
odir = os.path.join(exp, 'RGI62', 'b_080', 'L5', 'summary')
title = exp.replace(idir, '').replace('/match_geod', '').replace('/elev_bands/', '/').replace('/', ' ')
for rgi_reg in np.arange(1, 20):
rgi_reg = '{:02d}'.format(rgi_reg)
f = odir + '/historical_run_output_extended_{}.nc'.format(rgi_reg)
try:
with xr.open_dataset(f) as ds:
ds = ds.sel(time=slice(2000, None))
area = ds.area.isel(time=0)
area = area.where(area > 0)
dmdtda = (ds.volume.isel(time=-1) - ds.volume.isel(time=0)) / area / len(ds.time) * 900 # ice density
except FileNotFoundError:
pass
dfz.loc[rgi_reg, 'SMB_OGGM_run'] = np.average(dmdtda.dropna(dim='rgi_id'), weights=area.dropna(dim='rgi_id'))
f, ax = plt.subplots()
dfz.plot(ax=ax, y='SMB_OGGM_run', marker='o', linestyle='none', markersize=9);
dfz.plot(ax=ax, y='SMB_HUGUONNET', yerr='SMB_HUGUONNET_err', marker='o', linestyle='none');
plt.xticks(range(0, len(dfz.index)), labels=dfz.index);
plt.title(title);