Glaciers explorer using Datashader#


This notebook provides an annotated HoloViews+Panel implementation of a dashboard originally developed by Fabien Maussion in Plotly+Dash for viewing data about the Earth’s glaciers from the Open Global Glacier Model. To run it, first:

conda install -c pyviz pandas=0.24 param=1.10.0 panel=0.10.1 holoviews=1.13.5 datashader=0.11.1

Next, save the data file as data/oggm_glacier_explorer.csv (and gzip it if desired).

The dashboard can then be used here as a cell in the Jupyter notebook, or you can run it as a separate server using:

panel serve glaciers.ipynb --show

This notebook is essentially the same as Glaciers.ipynb but uses unaggregated data that is practical only with Datashader.

import numpy as np
import pandas as pd
import holoviews as hv
import datashader as ds
import panel as pn

from colorcet import bmy
from holoviews.operation.datashader import rasterize, datashade


Load the data#

Here we will load the glaciers data and project the latitudes and longitudes to Google Mercator coordinates, which will allow us to plot it on top of a tile source. We use the pn.state.as_cached function to cache the data to ensure that only the first visitor to our app has to load the data.

def load_data():
    df = pd.read_csv('data/oggm_glacier_explorer.csv')
    df['latdeg'] = df.cenlat
    df['x'], df['y'] = ds.utils.lnglat_to_meters(df.cenlon, df.cenlat)
    return df

df = pn.state.as_cached('glaciers', load_data)

rgi_id cenlon cenlat area_km2 glacier_type terminus_type mean_elev max_elev min_elev avg_temp_at_mean_elev avg_prcp latdeg x y
213745 RGI60-18.03533 170.354 -43.4215 0.189 Glacier Land-terminating 1704.858276 2102.0 1231.0 2.992555 6277.991881 -43.4215 1.896372e+07 -5.376350e+06
213746 RGI60-18.03534 170.349 -43.4550 0.040 Glacier Land-terminating 2105.564209 2261.0 1906.0 0.502311 6274.274146 -43.4550 1.896316e+07 -5.381486e+06
213747 RGI60-18.03535 170.351 -43.4400 0.184 Glacier Land-terminating 1999.645874 2270.0 1693.0 1.187901 6274.274146 -43.4400 1.896339e+07 -5.379186e+06
213748 RGI60-18.03536 170.364 -43.4106 0.111 Glacier Land-terminating 1812.489014 1943.0 1597.0 2.392771 6154.064456 -43.4106 1.896483e+07 -5.374680e+06
213749 RGI60-18.03537 170.323 -43.3829 0.085 Glacier Land-terminating 1887.771484 1991.0 1785.0 1.351039 6890.991816 -43.3829 1.896027e+07 -5.370436e+06

Plot the data#

As you can see in the dataframe, there are a lot of things that could be plotted about this dataset, but following the previous version let’s focus on the lat/lon location, elevation, temperature, and precipitation. We’ll use tools from HoloViz, starting with HoloViews as an easy way to build interactive Bokeh plots. So that we can use the full glacier database with good performance, we’ll have Datashader pre-render some of the plots as images before they reach the browser.

To start, let’s declare a HoloViews object that captures English-text descriptions of the various columns in the dataframe, in a way that subsequent plots can all inherit without having to repeat that information:

data = hv.Dataset(df, [('x', 'Longitude'), ('y', 'Latitude')],
                     [('avg_prcp', 'Annual Precipitation (mm/yr)'),
                      ('area_km2', 'Area'), ('latdeg', 'Latitude (deg)'),
                      ('avg_temp_at_mean_elev', 'Annual Temperature at avg. altitude'), 
                      ('mean_elev', 'Elevation')])
total_area = df.area_km2.sum()

print(data, len(data), total_area)
:Dataset   [x,y]   (avg_prcp,area_km2,latdeg,avg_temp_at_mean_elev,mean_elev) 213750 613225.6620000001

Here we’ve declared that x and y (the projected lat,lon location of the center of the glacier) are the “key dimensions” (independent values that specify which glacier this is), and the rest are “value dimensions” (various dependent values characterizing that particular sample).

Next, let’s define various options that will control the appearance of our plots:

geo_kw    = dict(aggregator=ds.sum('area_km2'), x_sampling=1000, y_sampling=1000)
elev_kw   = dict(cmap='#7d3c98')
temp_kw   = dict(num_bins=50, adjoin=False, normed=False, bin_range=data.range('avg_temp_at_mean_elev'))
prcp_kw   = dict(num_bins=50, adjoin=False, normed=False, bin_range=data.range('avg_prcp'))

size_opts = dict(min_height=400, min_width=600, responsive=True)
geo_opts  = dict(size_opts, cmap=bmy, logz=True, colorbar=True, xlabel='', ylabel='')
elev_opts = dict(size_opts, show_grid=True)
temp_opts = dict(size_opts, fill_color='#f1948a', default_tools=[], toolbar=None, alpha=1.0)
prcp_opts = dict(size_opts, fill_color='#85c1e9', default_tools=[], toolbar=None, alpha=1.0)

Using these options with HoloViews, we can plot various combinations of the variables of interest:

geo_bg = hv.element.tiles.EsriImagery().opts(alpha=0.6, bgcolor="black")
geopoints = hv.Points(data, vdims=['area_km2']).opts(**geo_opts)

(geo_bg*rasterize(geopoints, **geo_kw).options(**geo_opts) + 
 datashade(, 'mean_elev','latdeg', []), **elev_kw).options(**elev_opts) + 
 data.hist('avg_temp_at_mean_elev', **temp_kw).options(**temp_opts) +
 data.hist('avg_prcp',              **prcp_kw).options(**prcp_opts)).cols(2)