Visualization/bokeh_holoviews_datashader.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive Visualization with Bokeh, HoloViews, and Datashader\n",
    "\n",
    "<br>Owner: **Keith Bechtol** ([@bechtol](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@bechtol))\n",
    "<br>Last Verified to Run: **2018-09-14**\n",
    "<br>Verified Stack Release: **v16.0, w201831, w201841**\n",
    "\n",
    "This notebook demonstrates a few of the interactive features of the [Bokeh](https://bokeh.pydata.org/en/latest/), [HoloViews](http://holoviews.org/), and [Datashader](http://datashader.org/) plotting packages in the notebook environment. These packages are part of the [PyViz](http://pyviz.org/) set of python tools intended for visualization use cases in a web browser, and can be used to create quite sophisticated dashboard-like interactive displays and widgets. The goal of this notebook is to provide an introduction and starting point from which to create more advanced, custom interactive visualizations. To get inspired, check out this beautiful [example notebook](https://github.com/lsst-sqre/notebook-demo/blob/master/experiments/QA-notebooks/coadd_9615_HSC-R.ipynb) using HSC data created with the [qa_explorer](https://github.com/timothydmorton/qa_explorer) tools.\n",
    "\n",
    "### Learning Objectives\n",
    "After working through and studying this notebook you should be able to\n",
    "   1. Use `bokeh` to create interactive figures with brushing and linking between multiple plots\n",
    "   2. Use `holoviews` and `datashader` to create two-dimensional histograms with dynamic binning to efficiently explore large datasets   \n",
    "\n",
    "Other techniques that are demonstrated, but not empasized, in this notebook are\n",
    "   1. Use `parquet` to efficiently access large amounts of data\n",
    "\n",
    "### Logistics\n",
    "This notebook is intended to be runnable on `lsst-lspdev.ncsa.illinois.edu` from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.\n",
    "\n",
    "Note that occasionally the notebook may seem to stall, or the interactive features may seem disabled. If this happens, usually a restart of the kernel fixes the issue. You might also need to log out of the LSP and start a \"large\" instance of the JupyterLab environment. In some examples shown in this notebook, the order in which the cells are run is important for understanding the interactive features, so you may want to re-run the set of cells in a given section if you encounter unexpected behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "You can find the Stack version by using `eups list -s` on the terminal command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the Stack am I using?\n",
    "! echo $HOSTNAME\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.io.fits as pyfits\n",
    "\n",
    "import bokeh\n",
    "from bokeh.io import output_file, output_notebook, show\n",
    "from bokeh.layouts import gridplot\n",
    "from bokeh.models import ColumnDataSource, Range1d, HoverTool, Selection\n",
    "from bokeh.plotting import figure, output_file\n",
    "\n",
    "import holoviews as hv\n",
    "from holoviews import streams\n",
    "from holoviews.operation.datashader import datashade, dynspread, rasterize\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need this line to display bokeh plots inline in the notebook\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prelude: Data Sample\n",
    "\n",
    "The data in the following example comes from the Dark Energy Survey Data Release 1 (DES DR1). The input data for this example obtained with the M2 globular cluster database query in Appendix C of the [DES DR1 paper](https://arxiv.org/abs/1801.03181) from the [DES Data Release page](https://des.ncsa.illinois.edu/releases/dr1/dr1-access)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "infile = '/project/kbechtol/des/dr1/dr1_m2_dered_test.fits'\n",
    "reader = pyfits.open(infile)\n",
    "data = reader[1].data\n",
    "reader.close()\n",
    "\n",
    "data = data[data['MAG_AUTO_G_DERED'] < 26.]\n",
    "print(len(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Brushing and linking between scatter plots with Bokeh\n",
    "\n",
    "First, an example with brushing and linking between two panels showing different repsentations of the same dataset. A selection applied to either panel will highlight the selected points in the other panel.\n",
    "\n",
    "Based on http://bokeh.pydata.org/en/latest/docs/user_guide/interaction/linking.html#linked-brushing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra_target, dec_target = 323.36, -0.82\n",
    "\n",
    "mag = data['MAG_AUTO_G_DERED']\n",
    "color = data['MAG_AUTO_G_DERED'] - data['MAG_AUTO_R_DERED']\n",
    "\n",
    "# create a column data source for the plots to share\n",
    "source = ColumnDataSource(data=dict(x0=data['RA'] - ra_target,\n",
    "                                    y0=data['DEC'] - dec_target,\n",
    "                                    x1=color,\n",
    "                                    y1=mag,\n",
    "                                    ra=data['RA'],\n",
    "                                    dec=data['DEC'],\n",
    "                                    coadd_object_id=data['COADD_OBJECT_ID']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a custom hover tool on both panels\n",
    "hover_left = HoverTool(tooltips=[(\"(RA,DEC)\", \"(@ra, @dec)\"),\n",
    "                                 (\"(g-r,g)\", \"(@x1, @y1)\"),\n",
    "                                 (\"coadd_object_id\", \"@coadd_object_id\")])\n",
    "hover_right = HoverTool(tooltips=[(\"(RA,DEC)\", \"(@ra, @dec)\"),\n",
    "                                  (\"(g-r,g)\", \"(@x1, @y1)\"),\n",
    "                                  (\"coadd_object_id\", \"@coadd_object_id\")])\n",
    "TOOLS = \"box_zoom,box_select,lasso_select,reset,help\"\n",
    "TOOLS_LEFT = [hover_left, TOOLS]\n",
    "TOOLS_RIGHT = [hover_right, TOOLS]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new plot and add a renderer\n",
    "left = figure(tools=TOOLS_LEFT, plot_width=500, plot_height=500, output_backend=\"webgl\",\n",
    "              title='Spatial: Centered on (RA, Dec) = (%.2f, %.2f)'%(ra_target, dec_target))\n",
    "left.circle('x0', 'y0', hover_color='firebrick', source=source,\n",
    "            selection_fill_color='steelblue', selection_line_color='steelblue',\n",
    "            nonselection_fill_color='silver', nonselection_line_color='silver')\n",
    "left.x_range = Range1d(0.3, -0.3)\n",
    "left.y_range = Range1d(-0.3, 0.3)\n",
    "left.xaxis.axis_label = 'Delta RA'\n",
    "left.yaxis.axis_label = 'Delta DEC'\n",
    "\n",
    "# create another new plot and add a renderer\n",
    "right = figure(tools=TOOLS_RIGHT, plot_width=500, plot_height=500, output_backend=\"webgl\",\n",
    "               title='CMD')\n",
    "right.circle('x1', 'y1', hover_color='firebrick', source=source,\n",
    "             selection_fill_color='steelblue', selection_line_color='steelblue',\n",
    "             nonselection_fill_color='silver', nonselection_line_color='silver')\n",
    "right.x_range = Range1d(-0.5, 2.5)\n",
    "right.y_range = Range1d(26., 16.)\n",
    "right.xaxis.axis_label = 'g - r'\n",
    "right.yaxis.axis_label = 'g'\n",
    "\n",
    "p = gridplot([[left, right]])\n",
    "\n",
    "# The plots can be exported as html files with data embedded\n",
    "#output_file(\"bokeh_m2_example.html\", title=\"M2 Example\")\n",
    "\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the hover tool to see information about individual datapoints (e.g., the `coadd_object_id`). This information should appear automatically as you hover the mouse over the datapoints. Notice the data points highlighted in red on one panel with the hover tool are also highlighted on the other panel. \n",
    "\n",
    "Next, click on the selection box icon (with a \"+\" sign) or the selection lasso icon found in the upper right corner of the figure.  Use the selection box and selection lasso to make various selections in either panel by clicking and dragging on either panel. The selected data points will be displayed in the other panel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducing HoloViews Linked Streams\n",
    "\n",
    "If we want to do subsequent calculations with the set of selected points, we can use HoloViews [linked streams](http://holoviews.org/user_guide/Custom_Interactivity.html) for custom interactivity. The following visualization is a modification of this [example](http://holoviews.org/reference/streams/bokeh/Selection1D_points.html).\n",
    "\n",
    "For this visualization, as in the example above, use the selection box and selection lasso to datapoints on the left panel. The selected points should appear in the right panel.\n",
    "\n",
    "Finally, notice that as you change the selection on the left panel, the mean x- and y-values for selected datapoints are shown in the title of right panel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output size=150\n",
    "%%opts Points [tools=['box_select', 'lasso_select']]\n",
    "\n",
    "# Declare some points\n",
    "points = hv.Points((data['RA'] - ra_target, data['DEC'] - dec_target))\n",
    "\n",
    "# Declare points as source of selection stream\n",
    "selection = streams.Selection1D(source=points)\n",
    "\n",
    "# Write function that uses the selection indices to slice points and compute stats\n",
    "def selected_info(index):\n",
    "    selected = points.iloc[index]\n",
    "    if index:\n",
    "        label = 'Mean x, y: %.3f, %.3f' % tuple(selected.array().mean(axis=0))\n",
    "    else:\n",
    "        label = 'No selection'\n",
    "    return selected.relabel(label).options(color='red')\n",
    "\n",
    "# Combine points and DynamicMap\n",
    "# Notice the interesting syntax used here: the \"+\" sign makes side-by-side panels\n",
    "points + hv.DynamicMap(selected_info, streams=[selection])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we access the indices of the selected datapoints. We could use these indices to select a subset of full sample for further examination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(selection.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intermission: Rapid Data Access with Parquet\n",
    "\n",
    "For the next example, we want to use a much larger dataset. Let's open up some data from Gata Data Release 2 (Gaia DR2) with Parquet. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import pandas as pd\n",
    "import pyarrow.parquet as pq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "infiles = sorted(glob.glob('/project/shared/data/gaia_dr2_1am/*.parquet'))\n",
    "print('There are %i total files in the directory'%(len(infiles)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df_array = []\n",
    "for ii in range(0, 10):\n",
    "    print(infiles[ii])\n",
    "    columns = ['ra', 'dec', 'phot_g_mean_mag'] # 'phot_g_mean_mag', 'phot_bp_mean_mag', 'phot_rp_mean_mag']\n",
    "    df_array.append(pq.read_table(infiles[ii], columns=columns).to_pandas())\n",
    "df = pd.concat(df_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Dataframe contains %.2f M rows'%(len(df) / 1.e6))\n",
    "print(df.columns.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Visualizing Larger Datasets with Datashader\n",
    "\n",
    "The interactive features of Bokeh work well with datasets up to a few tens of thousands of data points. To efficiently explore larger datasets, we'd like to use another visualization model that offers better scalability, namely [Datashader](http://datashader.org/).\n",
    "\n",
    "In the examples below, notice that as one zooms in on the datashaded two-dimensional histograms, the bin sizes are dynamically adjusted to show finer or coarser granularity in the distribution. This allows one to interactively explore large datasets without having to manually adjust the bin sizes while panning and zooming. Zoom in all the way and you can see individual points (i.e., bins contain either zero or one count). If you zoom in far enough, the individual points are represented by extremely small pixels in datashader that are difficult to see. A solution is to `dynspread` instead of `datashade`, which will preserve a finite size of the plotted points.\n",
    "\n",
    "In this particular example, as we zoom in, we can see that the Gaia dataset has been sharded into narrow stripes in declination.\n",
    "\n",
    "The next cell also uses the concept of linked Streams in HoloViews for [custom interactivity](http://holoviews.org/user_guide/Custom_Interactivity.html), in this case to create a selection box. We'll use that selection box tool in the following cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output size=150 # This cell controls the size of the figure in the notebook\n",
    "#%%opts Points [tools=['box_select']]\n",
    "points = hv.Points((df.ra, df.dec)) # Create a holoviews object to hold and plot data\n",
    "#points = hv.Points(np.random.multivariate_normal((0, 0), [[1, 0.1], [0.1, 1]], (1000,))) # If you wanted a simple synthetic dataset\n",
    "\n",
    "# Create the linked streams instance\n",
    "boundsxy = (0, 0, 0, 0)\n",
    "box = streams.BoundsXY(source=points, bounds=boundsxy)\n",
    "bounds = hv.DynamicMap(lambda bounds: hv.Bounds(bounds), streams=[box]) \n",
    "\n",
    "# Apply the datashader\n",
    "dynspread(datashade(points, cmap=bokeh.palettes.Viridis256)) * bounds\n",
    "# The \"*\" syntax puts multiple plot elements on the same panel\n",
    "#datashade(points, cmap=bokeh.palettes.Viridis256) * bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we add callback functionality to the plot above and retrieve the indices of the selected points. First, use the box selection tool to create a selection box for the two-dimensional histogram above. Then run the cell below to count the number of datapoints within the selection region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selection = (points.data.x > box.bounds[0]) \\\n",
    "    & (points.data.y > box.bounds[1]) \\\n",
    "    & (points.data.x < box.bounds[2]) \\\n",
    "    & (points.data.y < box.bounds[3])\n",
    "print('The selection box contains %i datapoints'%(np.sum(selection)))\n",
    "if np.sum(selection) > 0:\n",
    "    print('\\nHere are some of the selected indices...')\n",
    "    print(np.nonzero(selection)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option is to make a second linked plot paired with the box selection on the two-dimensional histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, create a holoviews dataset instance. Here we label some of the columns.\n",
    "kdims = [('ra', 'RA(deg)'), ('dec', 'Dec(deg)')]\n",
    "vdims = [('phot_g_mean_mag', 'G(mag)')]\n",
    "ds = hv.Dataset(df, kdims, vdims)\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = hv.Points(ds)\n",
    "\n",
    "#boundsxy = (0, 0, 0, 0)\n",
    "boundsxy = (np.min(ds.data['ra']), np.min(ds.data['dec']), np.max(ds.data['ra']), np.max(ds.data['dec']))\n",
    "box = streams.BoundsXY(source=points, bounds=boundsxy)\n",
    "box_plot = hv.DynamicMap(lambda bounds: hv.Bounds(bounds), streams=[box])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function defines the custom callback functionality to update the linked histogram\n",
    "def update_histogram(bounds=bounds):\n",
    "    \n",
    "    selection = (ds.data['ra'] > bounds[0]) & \\\n",
    "                (ds.data['dec'] > bounds[1]) & \\\n",
    "                (ds.data['ra'] < bounds[2]) & \\\n",
    "                (ds.data['dec'] < bounds[3])\n",
    "    \n",
    "    selected_mag = ds.data.loc[selection]['phot_g_mean_mag']\n",
    "    \n",
    "    frequencies, edges = np.histogram(selected_mag)\n",
    "    \n",
    "    hist = hv.Histogram((np.log(frequencies), edges))\n",
    "    return hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output size=150\n",
    "dmap = hv.DynamicMap(update_histogram, streams=[box])\n",
    "datashade(points, cmap=bokeh.palettes.Viridis256) * box_plot + dmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that when you select different regions of the left panel with the box select tool, the histogram on the right is updated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Images\n",
    "\n",
    "The next example demonstrates image visualization at the pixel level with datashader."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First get a sensor image\n",
    "import lsst.daf.persistence as dafPersist\n",
    "\n",
    "butler = dafPersist.Butler('/project/shared/data/Twinkles_subset/output_data_v2')\n",
    "subset = butler.subset('calexp')\n",
    "dataid = subset.cache[0]\n",
    "image = butler.get('calexp', immediate=True, dataId=dataid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%opts Image  [height=600 width=650]\n",
    "%%opts Bounds (color='white')\n",
    "#%%output size=200\n",
    "\n",
    "# Use an actual sensor image\n",
    "bounds_img = (0, 0, image.getDimensions()[0], image.getDimensions()[1])\n",
    "img = hv.Image(np.log10(image.image.array), bounds=bounds_img).options(colorbar=True, cmap=bokeh.palettes.Viridis256) # Another option is logz=True\n",
    "\n",
    "boundsxy = (0, 0, 0, 0)\n",
    "box = streams.BoundsXY(source=img, bounds=boundsxy)\n",
    "bounds = hv.DynamicMap(lambda bounds: hv.Bounds(bounds), streams=[box])\n",
    "\n",
    "rasterize(img) * bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with the histograms, it is possible to use interactive callback features on the image plots, such as the selection box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's another version of the image with a tap stream instead of box select. Click on the image to place an 'X' marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%opts Image  [height=600 width=650]\n",
    "%%opts Points (color='white' marker='x' size=20)\n",
    "\n",
    "posxy = hv.streams.Tap(source=img, x=0.5 * image.getDimensions()[0], y=0.5 * image.getDimensions()[1])\n",
    "marker = hv.DynamicMap(lambda x, y: hv.Points([(x, y)]), streams=[posxy])\n",
    "\n",
    "rasterize(img) * marker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "'X' marks the spot! What's the value at that location? Execute the next cell to find out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The value at position (%.3f, %.3f) is %.3f'%(posxy.x, posxy.y, img[int(posxy.x), int(posxy.y)]))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}