opt = 'ssp'
opt = 'rcp'
#language = 'en'
language = 'de'
import os, numpy as np, pandas as pd
import bokeh
import holoviews as hv
from bokeh.models import HoverTool, CustomJSHover
import xarray as xr
from oggm.utils import tolist
from collections import OrderedDict
import panel as pn
hv.extension('bokeh')
pn.__version__ , hv.__version__, bokeh.__version__
('0.14.2', '1.15.3', '2.4.3')
if opt =='ssp':
from text import trads_ssp as trads
else:
from text import trads
from text import trads_temp_change
Data sources:
[January 10th 2023]. - data published for the following paper: Rounce, D. R., Hock, R., Maussion, F., Hugonnet, R., Kochtitzky, W., Huss, M., ... & McNabb, R. W. (2023). Global glacier change in the 21st century: Every increase in temperature matters. Science, 379(6627), 78-83. DOI: 10.1126/science.abo1324, https://www.science.org/doi/10.1126/science.abo1324
data = pd.read_hdf('./data/rgi62_era5_itmix_country_df.h5', 'df')
data['vol_itmix_km3'] = data['vol_itmix_m3'] * 1e-9
data
GLIMSId | BgnDate | EndDate | CenLon | CenLat | O1Region | O2Region | Area | Zmin | Zmax | ... | era5_avg_pcp | era5_elev | era5_landmask | era5_trend | era5_trend_std_err | era5_trend_p_values | vol_itmix_m3 | vol_bsl_itmix_m3 | Country_Name | vol_itmix_km3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
RGIId | |||||||||||||||||||||
RGI60-01.00001 | G213177E63689N | 20090703 | -9999999 | -146.8230 | 63.6890 | 01 | 02 | 0.360 | 1936 | 2725 | ... | 763.129312 | 1520.692749 | 0.993423 | 0.144383 | 0.131552 | 0.279312 | 7.638771e+06 | 0.00000 | United States of America | 0.007639 |
RGI60-01.00002 | G213332E63404N | 20090703 | -9999999 | -146.6680 | 63.4040 | 01 | 02 | 0.558 | 1713 | 2144 | ... | 801.464239 | 1472.621460 | 0.995941 | 0.159570 | 0.129708 | 0.226173 | 1.697646e+07 | 0.00000 | United States of America | 0.016976 |
RGI60-01.00003 | G213920E63376N | 20090703 | -9999999 | -146.0800 | 63.3760 | 01 | 02 | 1.685 | 1609 | 2182 | ... | 823.346929 | 1494.963135 | 0.995392 | 0.133943 | 0.130580 | 0.311494 | 5.969346e+07 | 0.00000 | United States of America | 0.059693 |
RGI60-01.00004 | G213880E63381N | 20090703 | -9999999 | -146.1200 | 63.3810 | 01 | 02 | 3.681 | 1273 | 2317 | ... | 823.346929 | 1494.963135 | 0.995392 | 0.133943 | 0.130580 | 0.311494 | 1.952248e+08 | 0.00000 | United States of America | 0.195225 |
RGI60-01.00005 | G212943E63551N | 20090703 | -9999999 | -147.0570 | 63.5510 | 01 | 02 | 2.573 | 1494 | 2317 | ... | 871.860850 | 1515.656616 | 0.998459 | 0.170513 | 0.132849 | 0.207087 | 1.221541e+08 | 0.00000 | United States of America | 0.122154 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
RGI60-19.02748 | G322268E53986S | 20020502 | -9999999 | -37.7325 | -53.9860 | 19 | 03 | 0.042 | 310 | 510 | ... | 1121.466184 | 6.950999 | 0.000000 | 0.086660 | 0.048094 | 0.079499 | 5.502906e+05 | 0.00000 | S. Geo. and the Is. | 0.000550 |
RGI60-19.02749 | G323864E54831S | 20030207 | -9999999 | -36.1361 | -54.8310 | 19 | 03 | 0.567 | 330 | 830 | ... | 1505.146630 | 236.044983 | 0.556025 | 0.129571 | 0.055983 | 0.026142 | 1.300672e+07 | 0.00000 | S. Geo. and the Is. | 0.013007 |
RGI60-19.02750 | G322698E54188S | 20030207 | -9999999 | -37.3018 | -54.1884 | 19 | 03 | 4.118 | 10 | 1110 | ... | 1264.378020 | 95.402214 | 0.358323 | 0.105411 | 0.053106 | 0.054408 | 2.506893e+08 | 359846.64917 | S. Geo. and the Is. | 0.250689 |
RGI60-19.02751 | G269573E68866S | 19870101 | -9999999 | -90.4266 | -68.8656 | 19 | 01 | 0.011 | 170 | 270 | ... | 667.099859 | 19.770145 | 0.046129 | -0.095884 | 0.113010 | 0.401494 | 1.068206e+05 | 0.00000 | United States of America | 0.000107 |
RGI60-19.02752 | G037714E46897S | 19660301 | -9999999 | 37.7140 | -46.8972 | 19 | 04 | 0.528 | 970 | 1170 | ... | 1205.984800 | 25.721743 | 0.067493 | 0.006869 | 0.042001 | 0.870965 | 1.489316e+07 | 0.00000 | South Africa | 0.014893 |
216502 rows × 41 columns
Select the Alps only:
data = data.loc[data.O1Region == '11']
data = data.loc[(data.CenLon > 2.5) & (data.CenLat > 43.2)]
data.plot(kind='scatter', x='CenLon', y='CenLat');
also read in the global temperature change in 2100 relative to 1850-1900 data:
pd_temp_ch = pd.read_csv('data/Rounce_et_al_2023/Global_mean_temp_deviation_2081_2100_rel_1850_1900.csv')
def temp_ch_ssp(ssp):
_temp_ch = np.round(pd_temp_ch.loc[pd_temp_ch.Scenario == ssp]['global_mean_deviation_degC'].mean(),1)
return f', ΔT~+{_temp_ch}°C'
_ds_med_vol_l = []
_ds_std_vol_l = []
if opt == 'rcp':
for rcp in ['rcp26', 'rcp45', 'rcp85']:
_ds = xr.open_dataset(f'data/Rounce_et_al_2023/R11_glac_mass_annual_50sets_2000_2100-{rcp}.nc')
# strangely NaN models have 0 values everywhere ...
gcms_real = _ds.model[_ds.sum(dim='glacier').sum(dim='year').glac_mass_annual != 0]
_ds = _ds.sel(model = gcms_real)
_ds = _ds.rename_dims({'glacier':'rgi_id'})
_ds = _ds.set_index({'rgi_id':'RGIId'})
print(f'amount of GCMs for {rcp} scenario: {len(_ds.model)}')
# I do here the mean although Rounce et al., uses the more robust median
_ds_med_vol = _ds.mean(dim='model').glac_mass_annual / 900 # kg -> m3
_ds_med_vol['rcp'] = rcp
_ds_med_vol_l.append(_ds_med_vol)
_ds_std_vol = _ds.std(dim='model').glac_mass_annual / 900 # kg -> m3
_ds_std_vol['rcp'] = rcp
_ds_std_vol_l.append(_ds_std_vol)
ds = xr.concat(_ds_med_vol_l, dim='rcp')
ds_uq = xr.concat(_ds_std_vol_l, dim='rcp')
ds.coords['rcp'] = ['RCP2.6', 'RCP4.5', 'RCP8.5']
ds_uq.coords['rcp'] = ['RCP2.6', 'RCP4.5', 'RCP8.5']
elif opt == 'ssp':
dict_nice_ssp = dict(zip(['ssp119','ssp126', 'ssp245', 'ssp370', 'ssp585'],
['SSP1-1.9', 'SSP1-2.6', 'SSP2-4.5', 'SSP3-7.0', 'SSP5-8.5']))
for ssp in ['ssp119','ssp126', 'ssp245', 'ssp370', 'ssp585']:
_ds = xr.open_dataset(f'data/Rounce_et_al_2023/R11_glac_mass_annual_50sets_2000_2100-{ssp}.nc')
# strangely NaN models have 0 values everywhere ...
gcms_real = _ds.model[_ds.sum(dim='glacier').sum(dim='year').glac_mass_annual != 0]
_ds = _ds.sel(model = gcms_real)
_ds = _ds.rename_dims({'glacier':'rgi_id'})
_ds = _ds.set_index({'rgi_id':'RGIId'})
_ds = _ds.dropna(dim='model')
print(f'amount of GCMs for {ssp} scenario: {len(_ds.model)}')
label_ssp = dict_nice_ssp[ssp] + temp_ch_ssp(ssp)
_ds_med_vol = _ds.mean(dim='model').glac_mass_annual / 900 # kg -> m3
_ds_med_vol['ssp'] = label_ssp
_ds_med_vol_l.append(_ds_med_vol)
_ds_std_vol = _ds.std(dim='model').glac_mass_annual / 900 # kg -> m3
_ds_std_vol['ssp'] = label_ssp
_ds_std_vol_l.append(_ds_std_vol)
ds = xr.concat(_ds_med_vol_l, dim='ssp')
ds_uq = xr.concat(_ds_std_vol_l, dim='ssp')
# we remove ssp119 as it only has 4 GCMs ...
ds = ds.isel(ssp=[1,2,3,4])
ds_uq = ds_uq.isel(ssp=[1,2,3,4])
#ds = ds.sel(rgi_id=data.index)
#ds_uq = ds_uq.sel(rgi_id=data.index)
print('\n')
print('these glaciers are not in Rounce et al., 2023, but in the RGI, as they melted away before 2022')
for r in data.index:
if r not in ds.rgi_id.values:
print(r)
data = data.drop(r)
amount of GCMs for rcp26 scenario: 10 amount of GCMs for rcp45 scenario: 10
/home/lilianschuster/anaconda3/envs/goggm_edu_html_apps/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /home/lilianschuster/anaconda3/envs/goggm_edu_html_apps/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
amount of GCMs for rcp85 scenario: 10 these glaciers are not in Rounce et al., 2023, but in the RGI, as they melted away before 2022
/home/lilianschuster/anaconda3/envs/goggm_edu_html_apps/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
RGI60-11.03235 RGI60-11.03237
# Mass sum for each of the RGI regions ...
_ds = xr.open_dataset('data/Rounce_et_al_2023/volume_rounce_reg_sum.nc').isel(ssp=[1,2,3,4])
_ds = _ds.rename_dims({'time':'year'})
_ds = _ds.set_index(year='time')
_ds.coords[opt] = ds[opt].values
ds_reg_m = (_ds.glac_mass_annual / 900).mean(dim='gcm')
ds_reg_std = (_ds.glac_mass_annual / 900).std(dim='gcm')
# common glaciers that are in "European Alps":
len(data)
3890
A multiline plot per selection of glaciers:
if opt == 'ssp':
ref_yr = 2022
hover = HoverTool(tooltips=[('Year', '@{year}'),
(f'Volume (% 2022)', '@{Volume (% 2022)}{int}'),
('+σ', '@{+σ}{int}'),
('-σ', '@{-σ}{int}'),
],
mode='mouse')
else:
ref_yr = 2017
hover = HoverTool(tooltips=[('Year', '@{year}'),
(f'Volume (% 2017)', '@{Volume (% 2017)}{int}'),
('+σ', '@{+σ}{int}'),
('-σ', '@{-σ}{int}'),
],
mode='mouse')
ds = ds.sel(year=slice(ref_yr,2100))
ds_uq = ds_uq.sel(year=slice(ref_yr,2100))
ds_reg_m = ds_reg_m.sel(year=slice(ref_yr,2100))
ds_reg_std = ds_reg_std.sel(year=slice(ref_yr,2100))
def sel_glaciers(rgi_ids):
if np.any(rgi_ids == 'Global'):
sel = ds_reg_m.sum(dim='rgi_region')
sel_t = (ds_reg_m + ds_reg_std).sum(dim='rgi_region')
sel_b = (ds_reg_m - ds_reg_std).sum(dim='rgi_region')
else:
sel = ds.sel(rgi_id=tolist(rgi_ids)).sum(dim='rgi_id')
sel_t = (ds + ds_uq).sel(rgi_id=tolist(rgi_ids)).sum(dim='rgi_id')
sel_b = (ds - ds_uq).sel(rgi_id=tolist(rgi_ids)).sum(dim='rgi_id')
sel_t = sel_t / sel.isel(year=0) * 100
sel_b = sel_b / sel.isel(year=0) * 100
sel = sel / sel.isel(year=0) * 100
return sel, sel_b, sel_t
def make_curve(rgi_ids, rcp=None, ssp = None,
line_dash=(0, 0), add_label=trads['sel_label'][language], color='#30a2da', opt='rcp',
short_label=False):
sel, sel_b, sel_t = sel_glaciers(rgi_ids)
if opt == 'rcp':
df = sel.sel(rcp=rcp).to_dataframe(name=f'Volume (% {ref_yr})')
df['+σ'] = sel_t.sel(rcp=rcp).to_series().clip(0)
df['-σ'] = sel_b.sel(rcp=rcp).to_series().clip(0)
label = add_label+rcp
elif opt=='ssp':
df = sel.sel(ssp=ssp).to_dataframe(name=f'Volume (% {ref_yr})')
df['+σ'] = sel_t.sel(ssp=ssp).to_series().clip(0)
df['-σ'] = sel_b.sel(ssp=ssp).to_series().clip(0)
if short_label:
label = add_label+ssp[:8]
else:
label = add_label+ssp
return hv.Curve(df, vdims=[f'Volume (% {ref_yr})', opt, '+σ', '-σ'],
kdims=['year'],
label=label).opts(tools=['hover'], line_dash=line_dash, color=color)
if opt == 'ssp':
# 'ssp119, ssp126', 'ssp245', 'ssp370', 'ssp585'
#colors_ssp = ['#00a9cf', '#003466',
# '#f69320', '#df0000', '#980002']
def sel_overlay(rgi_ids, line_dash=(0, 0), add_label=trads['sel_label'][language], short_label=False):
# use IPCC SSP colors
return hv.Overlay(
#make_curve(rgi_ids, ssp='SSP1-1.9', line_dash=line_dash, add_label=add_label, color='#00a9cf', opt='ssp') *
make_curve(rgi_ids, ssp='SSP1-2.6'+temp_ch_ssp('ssp126'),
line_dash=line_dash, add_label=add_label, color='#003466', opt='ssp', short_label=short_label) *
make_curve(rgi_ids, ssp='SSP2-4.5'+temp_ch_ssp('ssp245'),
line_dash=line_dash, add_label=add_label, color='#f69320', opt='ssp', short_label=short_label) *
make_curve(rgi_ids, ssp='SSP3-7.0'+temp_ch_ssp('ssp370'),
line_dash=line_dash, add_label=add_label, color='#df0000', opt='ssp', short_label=short_label) *
make_curve(rgi_ids, ssp='SSP5-8.5'+temp_ch_ssp('ssp585'),
line_dash=line_dash, add_label=add_label, color='#980002', opt='ssp', short_label=short_label)
)
else:
def sel_overlay(rgi_ids, line_dash=(0, 0), add_label=trads['sel_label'][language], short_label=np.NaN):
return hv.Overlay(
make_curve(rgi_ids, rcp='RCP2.6', line_dash=line_dash, add_label=add_label, color='#30a2da') *
make_curve(rgi_ids, rcp='RCP4.5', line_dash=line_dash, add_label=add_label, color='#e5ae38') *
make_curve(rgi_ids, rcp='RCP8.5', line_dash=line_dash, add_label=add_label, color='#fc4f30')
)
Define the regions:
curve_dict = OrderedDict()
curve_dict[trads['alps_label'][language]] = sel_overlay(data.index, short_label=True)
for c in data.Country_Name.unique():
if c == 'Slovenia':
continue
curve_dict[c] = sel_overlay(data.loc[data.Country_Name==c].index, short_label=True)
# get rgi ids of glaciers from Pitztal, Ötztal, and Stubaital (new since Oct 2023)
glac_piz_otz_stb = pd.read_csv('data/RGI_PIZ_OTZ_STB_created_by_Patrick_Schmitt.csv').rgi_id.values
curve_dict['Oetztal/Pitztal/Stubaital'] = sel_overlay(list(glac_piz_otz_stb))
# below is the WRONG old selection:
#curve_dict['Oetztal/Pitzal'] = sel_overlay(['RGI60-11.00670', 'RGI60-11.00666', 'RGI60-11.00663', 'RGI60-11.00648', 'RGI60-11.00674'], short_label=True)
#curve_dict['Global'] = sel_overlay('Global')
# This could be added at whish
# curve_dict['Brunnenkogelferner'] = sel_overlay('RGI60-11.00670')
# curve_dict['Mittelberferner'] = sel_overlay('RGI60-11.00666')
# curve_dict['Karlesferner'] = sel_overlay('RGI60-11.00663')
# curve_dict['Rettenbachferner'] = sel_overlay('RGI60-11.00648')
# curve_dict['Tiefenbachferner'] = sel_overlay('RGI60-11.00674')
hmap = hv.HoloMap(curve_dict, kdims='Region', sort=False)
fplot = sel_overlay(data.index, line_dash=(4, 4),
add_label=trads['alps_label'][language]+': ').opts(tools=['hover']) * hmap.opts(tools=['hover'])
fplot = fplot.opts(width=900, height=500, fontsize={'title': 16, 'labels': 14, 'xticks': 12, 'yticks': 12, 'legend':12},
ylim=(-0.02,106), #yaxis='right',
#legend_position='right',
#legend_cols=True)
)
# would be nice to have a secondary y axis on the right side ... secondary_yaxis="right")
# to have the legend outside or at least two columns and to have another
title = pn.pane.Markdown(sizing_mode='stretch_height', width=1080)
title.object = '<div style="font-size:38px; color: #326a86; font-weight: bold" >{}</div>'.format(trads['Title'][language])
#where to find this logo?
oggm_logo = '<a href="http://edu.oggm.org"><img src="https://raw.githubusercontent.com/zschirmeister/glacier-gallery/master/oggm_loupe.png" width=220></a>'
pn_logo = '<a href="https://panel.pyviz.org"><img src="https://panel.pyviz.org/_static/logo_stacked.png" width=46 height=39></a>'
holo_logo = '<a href="https://holoviz.org/"><img src="https://raw.githubusercontent.com/pyviz/holoviews/master/doc/_static/logo.png" width=46 height=39></a>'
intro_text = trads['instruction'][language]
instruction = pn.pane.Markdown('<span style="color:#326a86;font-size:15px">'+intro_text+'</span>')
#sel_region = '<div>{}</div>'.format(trads['sel_region'][language])
#logos = pn.Row(pn_logo, holo_logo)
# Rounce, D.R., R. Hock, and F. Maussion. 2022. Global PyGEM-OGGM Glacier Projections with RCP and SSP Scenarios, Version 1. [R11_glac_mass_annual]. Boulder, Colorado USA.
#NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/P8BN9VO9N5C7.
#[January 10th 2023].
header = pn.Column(pn.Pane(oggm_logo))
source = pn.pane.Markdown('Data source: [Rounce, D.R., R. Hock, and F. Maussion. 2022. Global PyGEM-OGGM Glacier Projections with RCP and SSP Scenarios](https://doi.org/10.5067/P8BN9VO9N5C7)',
width=500)
app = pn.Column(header, pn.Spacer(height=15), title, pn.Spacer(height=5), instruction, fplot, source)
# app.servable(title='Future glacier evolution in the Alps')
# transform app to model and save it as html-document:
app_to_save = pn.panel(app)
app_to_save.save(f'alps_future-app_rounce_{opt}_{language}.html', embed=True, title=trads['Title'][language])