In [ ]:
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()

Read the data

In [ ]:
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()
In [ ]:
# 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
In [ ]:
total_area = data.area.sum()
total_vol = data.vol_itmix_km3.sum()

Convert to HoloViews dataset

In [ ]:
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))

Sea-level equivalent computation

In [ ]:
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)

Plots kwargs

In [ ]:
# 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'))
In [ ]:
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)
In [ ]:
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)

Language settings

In [ ]:
from international import trads, supported_languages
language = 'en'

Plots

In [ ]:
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'])))])
In [ ]:
static_slr  = base_slr(data)
In [ ]:
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

In [ ]:
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)

Dynamic selection

In [ ]:
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)

Dynamic plots

In [ ]:
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
In [ ]:
# set range
elevation = elevation.redim(latdeg=dict(range=(-90, 90)), Zmed=dict(range=(0, 8000)))
In [ ]:
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()

Language selector

In [ ]:
def language_selector_function(arg=None):
    global language
    language = language_selector.value
    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=supported_languages,
                                             inline=True,
                                             margin=(0, 0),
                                             width=70)
language_selector.param.watch(language_selector_function, 'value');

functions to change language with datashader

In [ ]:
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]

Put everything together

In [ ]:
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)
In [ ]:
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)
              )
In [ ]:
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)
In [ ]:
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]))
In [ ]:
top = pn.Row(left, pn.layout.Spacer(width=40), overview, pn.layout.Spacer(width=30), geomap.options(hooks=[geo_language]))
In [ ]:
app = pn.Column(pn.Row(pn.layout.Spacer(width=1300), language_selector),
                top,
                pn.layout.Spacer(height=5),
                plots,
                explanation)
In [ ]:
app.servable(title='World glaciers explorer')
In [ ]: