08b_Interactive_Catalog_Visualization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250, style=\"padding: 10px\"> \n",
    "<p><p><p><p><p><p>\n",
    "<b>Interactive Catalog Visualization</b> <br>\n",
    "Last verified to run on 2022-05-31 with LSST Science Pipelines release w_2022_22<br>\n",
    "Contact authors: Leanne Guy <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "Minimum Container Size: medium <br>\n",
    "Questions welcome at <a href=\"https://community.lsst.org/c/support/dp0\">community.lsst.org/c/support/dp0</a> <br>\n",
    "Find DP0 documentation and resources at <a href=\"https://dp0-1.lsst.io\">dp0-1.lsst.io</a> <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial was inspired by a notebook originally developed by Keith Bechtol in the context of the LSST Stack Club. It has been updated and extended for DP0.1 by Leanne Guy. Please consider acknowledging Leanne Guy and Keith Bechtol in any publications or software releases that make use of this notebook's contents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "This tutorial, together with tutorial `08a_Interactive_Image_Visualization`, introduces three open-source Python libraries that enable powerful interactive visualization of catalogs.\n",
    " 1. [**HoloViews**](http://holoviews.org): Produce high-quality interactive visualizations easily by annotating plots and images rather than using direct calls to a plotting library\n",
    " 2. [**Bokeh**](https://bokeh.org): A powerful data visualization library that provides interactive tools including brushing and linking between multiple plots. `Holoviews` + `Bokeh`\n",
    " 3. [**Datashader**](https://datashader.org): Accurately render very large datasets quickly and flexibly.\n",
    " \n",
    "These packages are part of the [Holoviz](http://holoviz.org/) ecosystem of tools intended for visualization in a web browser and can be used to create quite sophisticated dashboard-like interactive displays and widgets. The goal of this tutorial is to provide an introduction and starting point from which to create more advanced, custom interactive visualizations. Holoviz has a [vibrant and active community](https://discourse.holoviz.org/) where you can ask questions and discuss vizualizations with a global community. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistics\n",
    "This notebook is intended to be runnable on `data.lsst.cloud`. 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 RSP 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. Note that some of the examples require manual selection of points on a graph to run correctly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# General python imports\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings\n",
    "\n",
    "# Update this option setting as you prefer\n",
    "pd.set_option('display.max_rows', 5)\n",
    "\n",
    "# Astropy\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.units.quantity import Quantity\n",
    "\n",
    "# LSST imports\n",
    "from lsst.rsp import get_tap_service\n",
    "\n",
    "# Bokeh and Holoviews for visualization\n",
    "import bokeh\n",
    "from bokeh.io import output_notebook, show\n",
    "from bokeh.models import ColumnDataSource, Range1d, HoverTool\n",
    "from bokeh.models import CDSView, GroupFilter\n",
    "from bokeh.plotting import figure, gridplot\n",
    "from bokeh.transform import factor_cmap\n",
    "\n",
    "import holoviews as hv\n",
    "from holoviews import streams, opts\n",
    "from holoviews.operation.datashader import datashade, dynspread\n",
    "from holoviews.plotting.util import process_cmap\n",
    "\n",
    "import datashader as dsh\n",
    "\n",
    "# Set the holoviews plotting library to be bokeh\n",
    "# You will see the holoviews + bokeh icons displayed when the library is loaded successfully\n",
    "hv.extension('bokeh')\n",
    "\n",
    "# Display bokeh plots inline in the notebook\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a message about \\\"Patching auth into notebook.base.handlers ...\\\" appeared above, it is ok to ignore, as are messages about the number of threads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What versions of bokeh and holoviews nd datashader are we working with?\n",
    "# This is important when referring to online documentation as\n",
    "# APIs can change between versions.\n",
    "print(\"Bokeh version: \" + bokeh.__version__)\n",
    "print(\"Holoviews version: \" + hv.__version__)\n",
    "print(\"Datashader version: \" + dsh.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ignore warnings\n",
    "import warnings\n",
    "from astropy.units import UnitsWarning\n",
    "warnings.simplefilter(\"ignore\", category=UnitsWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the LSST Science Pipelnes are we using?\n",
    "! echo $IMAGE_DESCRIPTION\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Data preparation\n",
    "\n",
    "The basis for any data visualization is the underlying data. In this tutorial we will work with tabular data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1 DP0.1 tabular dataset\n",
    "We will execute a cone search about a defined coordinate with a specified radius using the Rubin TAP service. For more details about using the TAP service and ADQL queries, please refer to tutorial `02_Intermediate_TAP_Query`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a Rubin TAP service instance. \n",
    "service = get_tap_service()\n",
    "assert service is not None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function to build a query passing a coordinate and a search radius\n",
    "def getQuery(c: SkyCoord, r: Quantity) -> str:\n",
    "    query = \"SELECT obj.ra, obj.dec, obj.objectId, obj.extendedness, \"\\\n",
    "            \"obj.mag_g_cModel, obj.mag_r_cModel, obj.mag_i_cModel, \"\\\n",
    "            \"truth_type, truth.match_objectId \" \\\n",
    "            \"FROM dp01_dc2_catalogs.object as obj \" \\\n",
    "            \"JOIN dp01_dc2_catalogs.truth_match as truth \" \\\n",
    "            \"ON truth.match_objectId = obj.objectId \" \\\n",
    "            \"WHERE CONTAINS(POINT('ICRS', obj.ra, obj.dec),\"\\\n",
    "            \"CIRCLE('ICRS', \" + str(c.ra.value) + \", \" + str(c.dec.value) + \", \" \\\n",
    "            + str(r.to(u.deg).value) + \" )) = 1 \" \\\n",
    "            \"AND obj.good = 1 \"  \\\n",
    "            \"AND truth.match_objectid >= 0 \" \\\n",
    "            \"AND truth.is_good_match = 1\"\n",
    "    return query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a reference position on the sky and a radius in degrees for a cone search\n",
    "c1 = SkyCoord(ra=55.65*u.degree, dec=-40.*u.degree, frame='icrs')\n",
    "radius = 4 * u.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = getQuery(c1, radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warnings about units that might appear above are ok to ignore."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This query will return a very large dataset, so we will use an asynchronous query. This could take a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "data = job.fetch_result().to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(data) == 23155236"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute some colours \n",
    "data['gmi'] = data['mag_g_cModel'] - data['mag_i_cModel']\n",
    "data['rmi'] = data['mag_r_cModel'] - data['mag_i_cModel']\n",
    "data['gmr'] = data['mag_g_cModel'] - data['mag_r_cModel']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map the truth type to a descriptive string\n",
    "# Catalog schema for the truth table can be found at:\n",
    "#  https://dp0-1.lsst.io/data-products-dp0-1/index.html#catalogs\n",
    "data['truth_type'] = data['truth_type'].map(\n",
    "    {1: 'galaxy', 2: 'star', 3: 'SNe'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert data[data[\"truth_type\"] == \"star\"].shape[0] == 504174\n",
    "assert data[data[\"truth_type\"] == \"galaxy\"].shape[0] ==  22651062\n",
    "assert data[data[\"truth_type\"] == \"SNe\"].shape[0] == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create a small subsample of this dataset of about 20000 points to demenstrate some of the basic functionality of the `holoviz` tools. In section 6, we will look at how to effectively visualize the full dataset of 200 million points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frac = 0.001\n",
    "data20K = data.sample(frac=frac, axis='index')\n",
    "assert len(data20K) == round(frac * len(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Holoviews\n",
    "\n",
    "[Holoviews](https://holoviews.org) supports easy analysis and visualization by annotating data rather than utilizing direct calls to plotting packages. For this tutorial, we will use [Bokeh](hrrps://bokeh.org) as the plotting library backend for Holoviews. This is defined in the `Setup` section above with the `hv.extension('bokeh')` call.  Holoviews supports several plotting libraries and there is an exercise to the user at the end of this section to explore using Holoviews with other plotting packages. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1 Visualizing tabular data with Holoviews\n",
    "\n",
    "The basic core primitives of Holoviews are [Elements](http://holoviews.org/Reference_Manual/holoviews.element.html); hv.Element. Elements are simple wrappers around your data that provide a semantically meaningful visual representation. An Element may be a set of Points, an Image, a Curve, a Histogram, etc. See the Holoviews [Reference Gallery](http://holoviews.org/reference/index.html) for all the various types of Elements that can be created with Holoviews. \n",
    "\n",
    "In this first example we will use the Holoviews [Scatter Element](http://holoviews.org/reference/elements/bokeh/Scatter.html) to quickly visualize the catalog data retrieved in section 1 as a scatter plot. HoloViews maintains a strict separation between content and presentation. This separation is achieved by maintaining sets of keyword values as `options` that specify how `Elements` are to appear.  In this first example we will apply the default options and remove the toolbar. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a simple scatter plot of the data using the Scatter element. \n",
    "hv.Scatter(data20K).options(toolbar=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `data20K` set contains many columns. If no columns are specified, explicitly the first 2 columns are taken for x and y respectively by the [Scatter Element](https://holoviews.org/reference/elements/bokeh/Scatter.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's bin the data in RA using the robust [Freedman Diaconis Estimator](https://numpy.org/doc/stable/reference/generated/numpy.histogram_bin_edges.html#numpy.histogram_bin_edges) and plot\n",
    "the resulting distribution using the Holoviews [Histogram Element](http://holoviews.org/reference/elements/bokeh/Histogram.html). \n",
    "We will also add in some basic plot options. Read more about about [customizing plots](https://holoviews.org/user_guide/Customizing_Plots.html) via `options`. Note that `options` can be shortened to `opts`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ra_bin, count) = np.histogram(data20K['ra'], bins='fd')\n",
    "ra_distribution = hv.Histogram((ra_bin, count)).opts(\n",
    "    title=\"RA distribution\",color='darkmagenta', \n",
    "    xlabel='RA', fontscale=1.2,\n",
    "    height=400, width=400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's create a layout of several plots. A `Layout` is a type of `Container` that can contain any HoloViews object. Other types of Containers that exist include `Overlay`, `Gridspace`, `Dynamicmap`, etc. See the Holoviews [Reference Gallery](http://holoviews.org/reference/index.html) for the full list of `Layouts` that can be created with Holoviews. See [Building Composite Objects](http://holoviews.org/user_guide/Building_Composite_Objects.html) for the full details about the ways Containers can be composed.\n",
    "\n",
    "The `+` operator is used to create a Layout. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slice the data and set some moore options\n",
    "skyplot = hv.Scatter(data20K[[\"ra\", \"dec\"]]).opts(\n",
    "    title=\"Skyplot\", toolbar='above', tools = ['hover'], \n",
    "    height=350, width=350, alpha=0.2, size=2) \n",
    "\n",
    "# Construct a layout using the `+` operator\n",
    "skyplots = skyplot + \\\n",
    "           ra_distribution.options(height=350, width=350)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skyplots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these two plots are not linked, they are two independent plots layed out next to each other. Try zooming in on the skyplot, you will notice that the data are not rebinned in the RA Distribution plot.  We will see how to link plots in Section 3. The tools however do apply to both plots. Try modifying both plots and then use the `reset` tool. Notice that both plots are reset to their original states. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's setup some default plot options to avoid duplicating long lists everytime we want to make a plot. As different plotting packages typically provide different customization capabilities, we will define one set of options for a Bokeh backend and one for a matplotlib backend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bokeh specific customizations as a python dictionary \n",
    "plot_style_bkh = dict(alpha=0.5, color='darkmagenta', \n",
    "                      marker='triangle', size=3,\n",
    "                      xticks=5, yticks=5,\n",
    "                      height=400, width=400, \n",
    "                      toolbar='above')\n",
    "\n",
    "# Matplotlib specific customizations\n",
    "plot_style_mpl = dict(alpha=0.2, color='c', marker='s', \n",
    "                      fig_size = 200, s=2, \n",
    "                      fontsize=14, xticks=8, yticks=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of subsetting a dataset to choose which columns to plot, Holoviews allows us to specifiy the dimensionality directly. \n",
    "`kdims` are the key dimensions or the independent variable(s) and `vdims` are the value dimensions or the dependenent variable(s). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the bokeh plot style\n",
    "plot_style = plot_style_bkh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Scatter(data20K, \n",
    "           kdims=['gmi'], vdims=['mag_g_cModel']\n",
    "          ).opts(invert_yaxis=True,\n",
    "                 xlabel=\" (g-i)\", ylabel=\"g\",\n",
    "                 **plot_style)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dimensions have be specified as strings above, but they are in fact rich objects. Dimension objects support a long descriptive label, which complements the short programmer-friendly name. Let's look at color-color diagram of the stars in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Axes as rich objects\n",
    "rmi = hv.Dimension('rmi', label='(r-i)', range = (-0.8,3.0))\n",
    "gmr = hv.Dimension('gmr', label='(g-r)', range = (-0.8,3.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a colour-colour scatter plot of just the stars in the dataset and also display the distribution of samples along both value dimensions using the `hist()` method of the [Scatter Element](http://holoviews.org/reference/elements/bokeh/Scatter.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stars = data20K[data20K[\"truth_type\"] == 'star']\n",
    "col_col = hv.Scatter(stars, kdims=gmr, \n",
    "                     vdims=rmi).opts(**plot_style)\n",
    "\n",
    "#  Use the hist method to show the distribution of samples along both value dimensions.\n",
    "col_col = col_col.hist(dimension = [gmr,rmi], \n",
    "                       num_bins=100, adjoin=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_col"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try zooming in on regions of the plot. The histograms are automatically recomputed.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The techniques to apply customizations in the cells above use standard Python syntax and are the recommended way to customize your visualizations in HoloViews. HoloViews also supports IPython magic commands. Magics are a much older approach that is not standard Python and is specific to notebooks. [HoloViews notebook magic](https://holoviews.org/user_guide/Notebook_Magics.html) supports both line and cell magics. Here is an example of using magics to plot the same spatial distribution of `Objects` as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%opts Scatter [tools=['hover'], toolbar='above',height=400, width=400](color='grey')\n",
    "hv.Scatter(data20K).opts(alpha=0.3, size=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our result set above contained a lot of columns. Often we want to be selective about which information we show in the hover tool and customize the names. We do this by creating a custom hover tool. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raDecHover = HoverTool(\n",
    "    tooltips=[\n",
    "        ( 'ra,dec', '@ra / @dec'),\n",
    "        ( 'rmag', '@mag_r_cModel'),\n",
    "        ( 'type', '@truth_type'),\n",
    "    ],\n",
    "    formatters={\n",
    "        'ra/dec' : 'printf',\n",
    "        'rmag' : 'numeral',\n",
    "        'type' : 'printf',\n",
    "    },\n",
    "    point_policy=\"follow_mouse\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Scatter(data20K).opts(tools=[raDecHover], \n",
    "                      **plot_style_bkh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It might be necessary to zoom in until points are distinguishable in order to notice how the hover tool box contents list only three properties per point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.0 Linked plots and data brushing with Bokeh\n",
    "\n",
    "A very useful feature of `Bokeh` is the ability to add connected interactivity between plots that show different attributes of the same data. This is called `linking`. With linked plots, we can carry out data `brushing`, whereby we select and maniplulate data synchronously across multiple linked plots. For example, if we link a skyplot with a colour-magnitude diagram of the same dataset, we can interactively explore the relationship between the positions of objects in each plot.  \n",
    "\n",
    "We will now look at the Bokeh plotting library directly to demonstrate how to set up 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",
    "This section is based on [Bokeh linked brushing](http://bokeh.pydata.org/en/latest/docs/user_guide/interaction/linking.html#linked-brushing)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1 Data preparation\n",
    "The basis for any data visualization is the underlying data. Getting the data preparation phase right is key to creating powerful visualizations. \n",
    "Bokeh works with a ColumnDataSource (CDS).  A CDS is essentially a collection of sequences of data that have their own unique column name. We will create a CDS from the data returned by the query above and pass it directly to bokeh. The CDS is the core of bokeh plots. Bokeh automatically creates a CDS from data passed as python lists or numpy arrays.  CDS are useful as they allow data to be shared between multiple plots and renderers, enabling brushing and linking. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a column data source for the plots to share. \n",
    "col_data = dict(x0=data20K['ra'] - c1.ra.value,\n",
    "                y0=data20K['dec'] - c1.dec.value,\n",
    "                x1=data20K['gmr'],\n",
    "                y1=data20K['mag_g_cModel'],\n",
    "                ra=data20K['ra'], dec=data20K['dec'])\n",
    "source = ColumnDataSource(data=col_data)\n",
    "\n",
    "# Additional data can be added to the CDS after creation\n",
    "source.data['objectId'] = data20K['objectId']\n",
    "source.data['rmi'] = data20K['rmi']\n",
    "source.data['gmr'] = data20K['gmr']\n",
    "source.data['mag_r_cModel'] = data20K['mag_r_cModel']\n",
    "source.data['truth_type'] = data20K['truth_type']\n",
    "source.data['extendedness'] = data20K['extendedness']\n",
    "\n",
    "# Create a view on truth_type stars\n",
    "stars = CDSView(source=source,\n",
    "                filters=[GroupFilter(column_name='truth_type', group=\"star\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2 Colour-Magnitude Diagram linked to a Skyplot\n",
    "\n",
    "We will use bokeh to plot a color-magnitude (g vs. g-i) diagram making use of the cModel magnitudes and a skyplot and then link the two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a custom hover tool for each panel. \n",
    "hover_left = HoverTool(tooltips=[(\"ObjectId\", \"@objectId\"),\n",
    "                                 (\"(ra,dec)\", \"(@ra, @dec)\"),\n",
    "                                 (\"type\", \"@truth_type\")\n",
    "                                 ])\n",
    "hover_right = HoverTool(tooltips=[(\"ObjectId\", \"@objectId\"),\n",
    "                                  (\"(g-r,g)\", \"(@x1, @y1)\"),\n",
    "                                  (\"extendedness\", \"@extendedness\")\n",
    "                                  ])\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. We will look at stars only\n",
    "stars = CDSView(source=source,\n",
    "                filters=[GroupFilter(column_name='truth_type', group=\"star\")])\n",
    "\n",
    "left = figure(tools=tools_left, \n",
    "              plot_width=400, plot_height=400,\n",
    "              title='Spatial: Centered on (RA, Dec) = (%.2f, %.2f)' %\n",
    "              (c1.ra.value, c1.dec.value))\n",
    "left.circle('x0', 'y0', hover_color='firebrick', \n",
    "            size=3, alpha=0.7,\n",
    "            source=source, \n",
    "            view=stars # select only stars\n",
    "           )\n",
    "left.x_range = Range1d(5.5, -5.5)\n",
    "left.y_range = Range1d(-5.5, 5.5)\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=400, plot_height=400, \n",
    "               title='CMD')\n",
    "right.circle('x1', 'y1', hover_color='firebrick', \n",
    "             size=4, alpha=0.8,\n",
    "             source=source, \n",
    "             view=stars  # Select only stars\n",
    "            )\n",
    "right.x_range = Range1d(-1.5, 2.8)\n",
    "right.y_range = Range1d(32., 16.)\n",
    "right.xaxis.axis_label = '(g-r)'\n",
    "right.yaxis.axis_label = 'g'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the grid of plots\n",
    "p = gridplot([[left, right]])\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the hover tool to see information about individual datapoints (e.g., the `ObjectId`). 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": [
    "### 4.0 Further analysis with Holoviews Linked Streams\n",
    "\n",
    "If we want to do subsequent calculations with the set of selected points, we can use HoloViews linked streams for custom interactivity. The following visualization is a modification of this example. As for 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. 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.\n",
    "\n",
    "Based on [Holoviews Selection1D points](http://holoviews.org/reference/streams/bokeh/Selection1D_points.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declare some points\n",
    "points = hv.Points((data20K['ra'] - c1.ra.value, data20K['dec'] - c1.dec.value)\n",
    "                  ).options(tools=['box_select', 'lasso_select'])\n",
    "\n",
    "# Declare points as source of selection stream\n",
    "selection = streams.Selection1D(source=points)\n",
    "\n",
    "# Define a function that uses the selection indices to slice points and compute stats\n",
    "def selected_info(index) -> str:\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",
    "\n",
    "# Combine points and DynamicMap\n",
    "# Notice the 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(len(selection.index))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.0  Visualizing Larger Datasets with Datashader"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\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",
    "The next cell also uses the concept of linked Streams in HoloViews for custom interactivity, in this case to create a selection box. We'll use that selection box tool in the following cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1 Color-color plot \n",
    "\n",
    "Here we plot a color-colour diagram of the cModel magnitudes obtained from the query in 1. Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create color-color plot using bokeh\n",
    "plot_options = {'plot_height': 400, 'plot_width': 800,\n",
    "                'tools': ['pan','box_zoom','box_select', 'wheel_zoom','reset', 'help']}\n",
    "\n",
    "p = figure(title=\"Colour-Colour Diagram (cModel magnitudes)\",\n",
    "           x_axis_label=\"(g-r)\", y_axis_label=\"(r-i)\",\n",
    "           x_range=(-2.0, 3.0), y_range=(-2.0, 3.0),\n",
    "           **plot_options)\n",
    "p.circle(x='gmr', y='rmi', source=source,\n",
    "         size=2, alpha=0.2,\n",
    "         hover_color='firebrick',\n",
    "         legend_field=\"truth_type\",\n",
    "         color=factor_cmap('truth_type', 'Category10_3',\n",
    "                           ['galaxy', 'star', 'SNe']))\n",
    "\n",
    "# Add a custiomized hover tool\n",
    "hover = HoverTool(tooltips=[(\"objectId\", \"@objectId\"),\n",
    "                            (\"(ra,dec)\", \"(@ra, @dec)\"),\n",
    "                            (\"(g-r,r-i)\", \"(@gmr, @rmi)\"),\n",
    "                            (\"type\", \"@truth_type\")])\n",
    "p.add_tools(hover)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that even with a medium sized dataset of ~20K points, this plot suffers from overplotting.  A classic strategy is to specify transparency of the glyphs so we can better see sparse and dense areas. In the plot above we have `alpha=0.2`. This helps,  but washes out the detail in the sparser regions. An additional problem is that we cannot add too many glyphs to any plot. \n",
    "\n",
    "Holoviews + Datashader allow us to plot millions to billions of points to produce much more informative plots. DataShader rasterizes or aggregates datasets into regular grids that can then be further analysed or viewed as plots or images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a holoviews object to hold and plot data\n",
    "points = hv.Points((source.to_df()['gmr'], \n",
    "                    source.to_df()['rmi'])).opts(\n",
    "    tools=['box_select', 'lasso_select'])\n",
    "\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",
    "p = dynspread(datashade(points, cmap=\"Viridis\"))\n",
    "p = p.opts(width=800, height=300,\n",
    "    padding=0.05, show_grid=True,\n",
    "    xlim=(-8.0, 18.0), ylim=(-18.0, 8.0),\n",
    "    xlabel=\"(g-r)\", ylabel=\"(r-i)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Render the datashaded plot\n",
    "p * bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `datashade` plot of the same color-color diagram as above does not require any magic-number parameters such as size and alpha and automatically ensures that there is no saturation or overplotting.  Select the `wheel zoom` and adjust the image as you interact with the plot. Note how the shades of color of the data points change according to the local density."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, this dataset is really too small to use datashader on. Let's see how to visualize the full 200 million dataset returned by the query with datashader. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Points Element for the data\n",
    "points = hv.Points((data['gmr'], \n",
    "                    data['rmi'])).opts(\n",
    "    tools=['box_select', 'lasso_select'])\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",
    "p = dynspread(datashade(points, cmap=\"Viridis\"))\n",
    "p = p.opts(width=800, height=300,\n",
    "    padding=0.05, show_grid=True,\n",
    "    xlim=(-8.0, 18.0), ylim=(-18.0, 8.0),\n",
    "    xlabel=\"(g-r)\", ylabel=\"(r-i)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Render the datashaded plot\n",
    "p * bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2 Adding a callback function\n",
    "Next we will add callback functionality to the colour-colour diagam above to retrieve the indices of selected points. We use the box selection tool to create a selection box for a two-dimensional histogram and then count the number of datapoints within the selection region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> STOP - Select some data points from the plot above using the box select tool before proceeding"
   ]
  },
  {
   "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)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3 Linked plots with interactive selection\n",
    "\n",
    "Now we will plot a spatial distribution on the sky of all the data (in the left-hand plot below), and link it to a two-dimensional histogram of the data (right-hand plot below). Use a box selection in the spatial distribution to change which data are included in the histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, create a holoviews dataset instance. \n",
    "# Here we label some of the columns.\n",
    "kdims = [('ra', 'RA(deg)'), ('dec', 'dec(deg)')]\n",
    "vdims = [('mag_r_cModel', 'r(mag)')]\n",
    "ds = hv.Dataset(data, kdims, vdims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = hv.Points(ds)\n",
    "boundsxy = (np.min(ds.data['ra']), np.min(ds.data['dec']),\n",
    "            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": [
    "# Custom callback functionality to update the linked histogram\n",
    "def log_inf(x) -> float:\n",
    "    return np.log(x) if x > 0 else 0\n",
    "\n",
    "\n",
    "def update_histogram(bounds=bounds) -> hv.Histogram:\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]['mag_r_cModel']\n",
    "    frequencies, edges = np.histogram(selected_mag)\n",
    "    hist = hv.Histogram((list(map(log_inf, frequencies)), edges))\n",
    "    return hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmap = hv.DynamicMap(update_histogram, streams=[box]).options(\n",
    "    height=400, width=400)\n",
    "datashade(points,\n",
    "          cmap=process_cmap(\"Viridis\", provider=\"bokeh\")) * \\\n",
    "box_plot.options(height=400, width=400) + \\\n",
    "dmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try changing the box selection and watch as the histogram is recomputed and displayed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally, delete the job\n",
    "job.delete()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.0  Optional exercises to the user \n",
    "\n",
    " 1. Holoviews works with a wide range of plotting libraries, Bokeh, matplotlib, plotly, mpld3, pygal to name a few. As an exercise, try changing the Holoviews plotting library to be `matplotlib` instead of `bokeh` in the `Setup` cell at the beginning of the notebook with `hv.extension('matplotlib')`. You will see the holoviews + matplotlib icons displayed when the library is loaded successfully. Run the cells in section 2.1 again and compare the outputs. Try again with some other plotting library. Don't forget to set the plotting library back to whichever you prefer to use for the rest of this tutorial.\n",
    "  \n",
    " 2. Try writing a different callback function to use in section 5.2. "
   ]
  }
 ],
 "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}