import os, numpy as np, pandas as pd, cartopy.crs as ccrs, bokeh
import holoviews as hv, geoviews as gv, datashader as ds, panel as pn
from holoviews.operation.datashader import rasterize, datashade
from colorcet import bmy, bgyw, isolum
from holoviews.util import Dynamic
from bokeh.models import HoverTool, CustomJSHover
hv.extension('bokeh', width=100)
pn.extension()
data = pd.read_hdf('./data/rgi62_era5_itmix_df.h5', 'df')
data = data.rename(columns={'Area':'area'})
data['latdeg'] = data.CenLat
data['vol_asl_m3'] = data.vol_itmix_m3 - data.vol_bsl_itmix_m3
data['vol_itmix_km3'] = data.vol_itmix_m3 * 1e-9
data.tail()
# Correct glaciers bsl
to_corr = data.loc[data.Zmed < 0]
data.loc[data.Zmed < 0, 'era5_avg_temp_at_zmed'] = to_corr['era5_avg_temp_at_zmed'] + to_corr.Zmed * 0.0065
data.loc[data.Zmed < 0, 'Zmed'] = 0
total_area = data.area.sum()
total_vol = data.vol_itmix_km3.sum()
data = gv.Points(data, [('CenLon', 'Longitude'), ('CenLat', 'Latitude')],
[('era5_avg_pcp', 'Annual Precipitation (mm/yr)'),
('area', 'Area'), ('latdeg', 'Latitude (deg)'),
('era5_avg_temp_at_zmed', 'Annual Temperature at avg. altitude (°C)'),
('Zmed', 'Mean elevation of the glacier (m a.s.l.)'), ('vol_asl_m3', 'Volume asl'),
('vol_bsl_itmix_m3', 'Volume bsl'),
('vol_itmix_km3', 'Volume km3'),
('era5_trend', 'Temperature trend 1979-2018 (°C per decade)')])
data = gv.Dataset(gv.operation.project_points(data))
def compute_slr(ice_vol_m3):
"""ice_vol in m³ gives slr in mm"""
rho = 900
A_oc = 362.5 * 1e12
return ice_vol_m3 * rho / A_oc
total_slr = compute_slr(data['vol_asl_m3']).sum()
np.testing.assert_allclose(total_slr, 332, atol=2)
# Datashader map
geo_kw = dict(aggregator=ds.sum('area'), x_sampling=1000, y_sampling=1000)
# Elev vs Lat scatter
elev_kw = dict(cmap='#7d3c98')
# Histograms
temp_kw = dict(num_bins=50, adjoin=False, normed=False, bin_range=data.range('era5_avg_temp_at_zmed'))
prcp_kw = dict(num_bins=50, adjoin=False, normed=False, bin_range=(0, 6000)) # for precip we crop large values
tren_kw = dict(num_bins=50, adjoin=False, normed=False, bin_range=data.range('era5_trend'))
size_opts_map = dict(height=520, width=715)
size_opts_his = dict(height=200, width=350)
size_opts_bar = dict(height=45, width=250)
size_opts_slr = dict(height=300, width=120)
size_text_bar = dict(height=20)
geo_opts = dict(size_opts_map, cmap=bmy, global_extent=False, logz=True, colorbar=True, colorbar_opts={'title':'Area (km²)'}, toolbar='above', projection=ccrs.GOOGLE_MERCATOR)
elev_opts = dict(size_opts_his, show_grid=True)
temp_opts = dict(size_opts_his, fill_color='#f1948a', default_tools=[], toolbar=None, ylabel='', alpha=1.0)
prcp_opts = dict(size_opts_his, fill_color='#85c1e9', default_tools=[], toolbar=None, ylabel='', alpha=1.0)
tren_opts = dict(size_opts_his, fill_color='#f4d34e', default_tools=[], toolbar=None, alpha=1.0)
slr_opts = dict(size_opts_slr, color='orange', default_tools=[], toolbar=None, xlabel='', show_legend=False, yticks=[0, 50, 100, 150,200,250,300,350,390], shared_axes=False)
glno_opts = dict(size_opts_bar, color='#326a86', default_tools=[], toolbar=None, alpha=0.8, invert_axes=True, show_legend=False, xaxis=None, yaxis=None, shared_axes=False)
area_opts = dict(size_opts_bar, color='#326a86', default_tools=[], toolbar=None, alpha=0.8, invert_axes=True, show_legend=False, xaxis=None, yaxis=None, shared_axes=False)
vol_opts = dict(size_opts_bar, color='#326a86', default_tools=[], toolbar=None, alpha=0.8, invert_axes=True, show_legend=False, xaxis=None, yaxis=None, shared_axes=False)
from international import trads, supported_languages
language = 'en'
def base_slr(data):
return hv.Bars([('asl', compute_slr(np.sum(data['vol_asl_m3']))), ('bsl', compute_slr(np.sum(data['vol_bsl_itmix_m3'])))])
static_slr = base_slr(data)
gl_number = len(data['area'])
def geo(data):
return gv.Points(data, crs=ccrs.PlateCarree).options(alpha=1)
def elev(data):
return data.to(hv.Scatter, 'Zmed', 'latdeg', [])
def temp(data):
return data.hist('era5_avg_temp_at_zmed', **temp_kw).options(**temp_opts)
def prcp(data):
return data.hist('era5_avg_pcp', **prcp_kw).options(**prcp_opts)
def tren(data):
return data.hist('era5_trend', **tren_kw).options(**tren_opts)
def slr(data):
slr_opts['ylabel'] = trads['bar_sealevel_y'][language]
return static_slr.opts(**slr_opts, alpha=0.1) * base_slr(data).opts(**slr_opts)
def gl_no(data):
return hv.Bars(('', len(data))).opts(**glno_opts)
def area(data):
return hv.Bars(('', np.sum(data['area']))).opts(**area_opts)
def vol(data):
return hv.Bars(('', np.sum(data['vol_itmix_km3']))).opts(**vol_opts)
def count1(data):
legend = trads['bar_glaciers_selected'][language]
text = '<p style="font-size:15px">{}</font>'
v1, v2 = len(data), gl_number
return hv.Div(text.format(legend).format(v1, v2)).options(**size_text_bar)
def count2(data):
legend = trads['bar_area'][language]
text = '<p style="font-size:15px">{}: {:.0f} km² ({:.1f}%)</font>'
v1, v2 = np.sum(data['area']), np.sum(data['area']) / total_area * 100
return hv.Div(text.format(legend, v1, v2)).options(**size_text_bar)
def count3(data):
legend = trads['bar_volume'][language]
text = '<p style="font-size:15px">{}: {:.0f} km³ ({:.1f}%)</font>'
v1, v2 = np.sum(data['vol_itmix_km3']), np.sum(data['vol_itmix_km3']) / total_vol * 100
return hv.Div(text.format(legend, v1, v2)).options(**size_text_bar)
def slr_text(data):
legend = trads['bar_sealevel_text'][language]
text = '<p style="font-size:15px">{}{:.1f} mm ({:.1f}%)</font>'
v1, v2 = compute_slr(np.sum(data['vol_asl_m3'])), compute_slr(np.sum(data['vol_asl_m3'])) / total_slr * 100
return hv.Div(text.format(legend, v1, v2)).options(**size_text_bar)
static_geo = rasterize(geo(data), **geo_kw).options(alpha=0.1, tools=['hover', 'box_select'], active_tools=['box_select'], **geo_opts)
static_elev = datashade(elev(data), **elev_kw).options(alpha=0.1, tools=[ 'box_select'], active_tools=['box_select'], toolbar=None, **elev_opts)
static_gl_no= gl_no(data).options(alpha=0.1)
static_area = area(data).options(alpha=0.1)
static_vol = vol(data).options(alpha=0.1)
static_tren = tren(data).options(alpha=0.1)
static_temp = temp(data).options(alpha=0.1)
static_prcp = prcp(data).options(alpha=0.1)
def combine_selections(**kwargs):
"""
Combines selections on all available plots into a single selection by index.
"""
if all(not v for v in kwargs.values()):
return slice(None)
selection = {}
for key, bounds in kwargs.items():
if bounds is None:
continue
elif len(bounds) == 2:
selection[key] = bounds
else:
xbound, ybound = key.split('__')
selection[xbound] = bounds[0], bounds[2]
selection[ybound] = bounds[1], bounds[3]
return sorted(set(data.select(**selection).data.index))
def select_data(**kwargs):
return data.iloc[combine_selections(**kwargs)] if kwargs else data
def clear_selections(arg=None):
geo_bounds.update(bounds=None)
elev_bounds.update(bounds=None)
temp_bounds.update(boundsx=None)
prcp_bounds.update(boundsx=None)
tren_bounds.update(boundsx=None)
Stream.trigger(selections)
from holoviews.streams import Stream, BoundsXY, BoundsX
geo_bounds = BoundsXY(source=static_geo, rename= {'bounds': 'CenLon__CenLat'})
elev_bounds = BoundsXY(source=static_elev, rename= {'bounds': 'Zmed__latdeg'})
temp_bounds = BoundsX( source=static_temp, rename= {'boundsx': 'era5_avg_temp_at_zmed'})
prcp_bounds = BoundsX( source=static_prcp, rename= {'boundsx': 'era5_avg_pcp'})
tren_bounds = BoundsX( source= static_tren, rename= {'boundsx': 'era5_trend'})
selections = [geo_bounds, elev_bounds, temp_bounds, prcp_bounds, tren_bounds]
dyn_data = hv.DynamicMap(select_data, streams=selections)
dyn_geo = rasterize(dyn_data.apply(geo), **geo_kw).options( **geo_opts)
dyn_elev = datashade(dyn_data.apply(elev), **elev_kw).options(**elev_opts)
dyn_temp = dyn_data.apply(temp)
dyn_prcp = dyn_data.apply(prcp)
dyn_count1= dyn_data.apply(count1)
dyn_count2= dyn_data.apply(count2)
dyn_count3= dyn_data.apply(count3)
dyn_slr_text = dyn_data.apply(slr_text)
dyn_tren = dyn_data.apply(tren)
dyn_slr = dyn_data.apply(slr)
dyn_gl_no = dyn_data.apply(gl_no)
dyn_area = dyn_data.apply(area)
dyn_vol = dyn_data.apply(vol)
geo_bg = gv.tile_sources.EsriImagery.options(alpha=1.0, bgcolor="white")
geomap = geo_bg * static_geo * dyn_geo
elevation = static_elev * dyn_elev
temperature = static_temp * dyn_temp
precipitation = static_prcp * dyn_prcp
trends = static_tren * dyn_tren
sealevelrise = dyn_slr
gl_num = static_gl_no * dyn_gl_no
area_bar = static_area * dyn_area
vol_bar = static_vol * dyn_vol
# set range
elevation = elevation.redim(latdeg=dict(range=(-90, 90)), Zmed=dict(range=(0, 8000)))
def set_button_name_language():
clear_button.name = trads['clear_button'][language]
clear_button = pn.widgets.Button(name='', width=170)
clear_button.param.watch(clear_selections, 'clicks')
set_button_name_language()
def language_selector_function(arg=None):
global language
language = language_selector.value.replace('汉语', 'cn')
change_language()
def change_hist_label_language():
temp_opts['xlabel'] = trads['temp_plot_x'][language]
prcp_opts['xlabel'] = trads['precip_plot_x'][language]
tren_opts['xlabel'] = trads['trend_plot_x'][language]
tren_opts['ylabel'] = trads['trend_plot_y'][language]
def change_language(arg=None):
set_instr_text()
set_explanation_text()
change_hist_label_language()
clear_selections()
set_button_name_language()
language_selector = pn.widgets.RadioBoxGroup(name='Select your language',
options=[l.replace('cn', '汉语') for l in supported_languages],
inline=True,
margin=(0, 0),
width=70)
language_selector.param.watch(language_selector_function, 'value');
def elev_language(plot, element):
plot.handles['xaxis'].axis_label = trads['elev_plot_x'][language]
plot.handles['yaxis'].axis_label = trads['elev_plot_y'][language]
def geo_language(plot, element):
plot.handles['xaxis'].axis_label = trads['map_plot_x'][language]
plot.handles['yaxis'].axis_label = trads['map_plot_y'][language]
oggm_logo = '<a href="https://edu.oggm.org"><img src="https://raw.githubusercontent.com/OGGM/world-glacier-explorer/master/img/logo_edu.png" width=180 height=79></a>'
fk_logo = '<a href="https://www.uibk.ac.at/foerderkreis1669/"><img src="https://raw.githubusercontent.com/OGGM/world-glacier-explorer/master/img/logo_1669.png" width=180 height=79></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>'
dasha_logo = '<a href="https://datashader.org/"><img src="https://raw.githubusercontent.com/pyviz/datashader/master/doc/_static/datashader-logo.png" width=46 height=39></a>'
logos = pn.Row(pn.layout.Spacer(width=10), pn_logo, holo_logo, dasha_logo)
left = pn.Column(oggm_logo, pn.layout.Spacer(height=1), logos, pn.layout.Spacer(height=1), fk_logo, pn.layout.Spacer(height=20), clear_button)
bars = pn.Row(pn.Column(pn.layout.Spacer(height=1),
pn.Pane(dyn_count1), pn.layout.Spacer(height=0), pn.Pane(gl_num, linked_axes=False),
pn.Pane(dyn_count2), pn.layout.Spacer(height=0), pn.Pane(area_bar, linked_axes=False),
pn.Pane(dyn_count3), pn.layout.Spacer(height=0), pn.Pane(vol_bar, linked_axes=False),
pn.Pane(dyn_slr_text), pn.layout.Spacer(height=1)),
pn.Pane(sealevelrise, linked_axes=False)
)
title = '<div style="font-size:35px">World glaciers explorer</div>'
instruction = pn.pane.Markdown(sizing_mode='stretch_width', height=100)
def set_instr_text():
instruction.object = trads['instructions'][language]
set_instr_text()
explanation = pn.pane.Markdown(sizing_mode='stretch_width', height=100)
def set_explanation_text():
explanation.object = trads['abbreviations'][language]
set_explanation_text()
overview = pn.Column(pn.Pane(title, width=400),
pn.Pane(instruction, width=470),
bars)
plots = pn.Row(pn.layout.Spacer(width=25), trends, pn.layout.Spacer(width=5), temperature, pn.layout.Spacer(width=5), precipitation, pn.layout.Spacer(width=5), elevation.options(hooks=[elev_language]))
top = pn.Row(left, pn.layout.Spacer(width=40), overview, pn.layout.Spacer(width=30), geomap.options(hooks=[geo_language]))
app = pn.Column(pn.Row(pn.layout.Spacer(width=1250), language_selector),
top,
pn.layout.Spacer(height=5),
plots,
explanation)
app.servable(title='World glaciers explorer')