08a_Interactive_Image_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 Image 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 `08b_Interactive_Catalog_Visualization`, introduces three open-source Python libraries that enable powerful interactive visualization of images and 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](https://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 warnings\n",
    "\n",
    "# Astropy\n",
    "from astropy.visualization import  ZScaleInterval, AsinhStretch\n",
    "\n",
    "# LSST imports\n",
    "from lsst.daf.butler import Butler\n",
    "\n",
    "# Bokeh and Holoviews for visualization\n",
    "import bokeh\n",
    "from bokeh.io import output_notebook\n",
    "from bokeh.models import HoverTool\n",
    "\n",
    "import holoviews as hv\n",
    "from holoviews import streams, opts\n",
    "from holoviews.operation.datashader import rasterize\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__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ignore warnings\n",
    "import warnings\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)"
   ]
  },
  {
   "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 images.  For DP0.1, images can only be accessed via the `Butler` (<a href=\"https://pipelines.lsst.io/modules/lsst.daf.butler/index.html\">documentation</a>), an LSST Science Pipelines software package that allows you to fetch the LSST data you want without having to know its location or format. For more details about using the Butler, please refer to tutorial `04_Intro_to_Butler`. For visualization examples with tabular data, see `08b_Interactive_Catalog_Visualization`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start by using the butler to retrieve a calexp image (specifying a visit, detector, and band) and a deep coadd image (specifying a tract, patch, and band)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate the Butler initializing it with the repository name and the DP0.1 collection identifier\n",
    "from lsst.daf.butler import Butler\n",
    "repo = 's3://butler-us-central1-dp01'\n",
    "collection = '2.2i/runs/DP0.1'\n",
    "butler = Butler(repo, collections=collection)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is ok to ignore a pink FutureWarning from the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a calibrated exposure and retrieve it via the Butler\n",
    "calexpId = {'visit': 192350, 'detector': 175, 'band': 'i'}\n",
    "calexp = butler.get('calexp', calexpId)\n",
    "assert calexp is not None\n",
    "print(f\"Visit: {calexp.info.getId()}, Band:{calexp.getFilterLabel().bandLabel}, Detector: {calexp.detector.getId()}\")\n",
    "\n",
    "# Note: This will trigger a warning from CFITSIO in w_2022_22. This warning can be safely ignored and will be corrected in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source table for this calexp exposure\n",
    "calexpSrc = butler.get('src', **calexpId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a deep coadded image and retrieve it via the Butler\n",
    "coaddId = {'tract': 4226, 'patch': 17, 'band': 'r'}\n",
    "coadd = butler.get('deepCoadd', coaddId)\n",
    "assert coadd is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get some infomation about the coadd, such as which visits went into constructing it. As an exercise, explore the information in the coaddInfo object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coaddInfo = coadd.getInfo()\n",
    "\n",
    "# Which visits went into constructing this coadd?\n",
    "coaddVisits = coaddInfo.getCoaddInputs().visits\n",
    "coaddVisits.asAstropy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source table for this coadd\n",
    "coaddSrc = butler.get('deepCoadd_forced_src', coaddId)"
   ]
  },
  {
   "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. \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, a Scatter, 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Visualizing exposure images with Holoviews\n",
    "\n",
    "In this first example we will use the Holoviews [Image Element](http://holoviews.org/reference/elements/bokeh/Image.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. \n",
    "\n",
    "In the tutorial `03_Image_Display_and_Manipulation` we saw how to use the `lsst.afw.display` library to visualize exposure images and in tutorial `03b_Image_Display_with_Firefly` we saw how to do the same using Firefly. In this example we demonstrate image visualization at the pixel level with Holoviews.\n",
    "\n",
    "We will use the holoviews Image Element to visualize a calexp. We will then overlay a Holoviews DynamicMap on the image to compute and display elements dynamically, allowing exploration of large datasets. DynamicMaps generate elements on the fly allowing exploration of parameters with arbitrary resolution. DynamicMaps are lazy in the sense they only compute as much data as the user wishes to explore. An Overlay is a collection of HoloViews objects that are displayed simultanously, e.g a Curve superimposed on a Scatter plot of data. You can build a Overlay between any two HoloViews objects, which can have different types using the * operator. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will use the `astropy.visualization` library to define an asinh stretch and zscale interval and apply them to the calexp object. These are the same transformations that were applied in `03_Image_Display_and_Manipulation` notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply a asinh/zscale mapping to the data \n",
    "transform = AsinhStretch() + ZScaleInterval()\n",
    "scaledImage = transform(calexp.image.array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LSST’s image classes (Image, Mask, MaskedImage, and Exposure) use a pixel indexing convention that is different from both the convention used by `numpy.ndarray` objects and the convention used in FITS images (as documented [here](https://pipelines.lsst.io/modules/lsst.afw.image/indexing-conventions.html)).  Most plotting tools assume pixel (0, 0) is in the upper left where we always assume (0,0) is in the lower left. Consequently, we flip the data array. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaledImage = np.flipud(scaledImage)\n",
    "bounds_img = (0, 0, calexp.getDimensions()[0], calexp.getDimensions()[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further details can be found at [Image Indexing, Array Views, and Bounding Boxes](https://pipelines.lsst.io/modules/lsst.afw.image/indexing-conventions.html) in the Rubin Science Pipelines and Data Products. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define some default plot options for the Image\n",
    "img_opts = dict(height=600, width=700, \n",
    "                xaxis=\"bottom\", \n",
    "                padding = 0.01, fontsize={'title': '8pt'},\n",
    "                colorbar=True, toolbar='right', show_grid=True,\n",
    "                tools=['hover']\n",
    "               )     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a function to autogenerate a plot title from the dataId.\n",
    "def dataIdToString(dataId: dict) -> str:\n",
    "    title = \"DC2 image: \"\n",
    "    for key, value in dataId.items():\n",
    "        title += str(key) + \": \" + str(value) + \" \"\n",
    "    return title.strip() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the Image element.\n",
    "img = hv.Image(scaledImage, bounds=bounds_img,\n",
    "               kdims=['x', 'y']).opts(\n",
    "    cmap = \"Greys_r\",  xlabel = 'X', ylabel ='Y',\n",
    "    title = dataIdToString(calexpId),\n",
    "    **img_opts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterize(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3 Overlaying source detections on an image\n",
    "\n",
    "Now let's overlay the source detections from the `Source` catalog on this image. We will use the Points Element for the detections to overlay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = calexpSrc.getX(), calexpSrc.getY()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f\"Number of src detections is, {len(coords[1])}\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Custom hover tool for the source detections\n",
    "detHoverTool = HoverTool(\n",
    "    tooltips=[\n",
    "        ( 'X', '@x{0.2f}'),\n",
    "        ( 'Y', '@y{0.2f}'),\n",
    "    ],\n",
    "    formatters={\n",
    "        'X' : 'printf',\n",
    "        'Y' : 'printf',\n",
    "    },\n",
    "    \n",
    ")\n",
    "detections = hv.Points(coords).opts(\n",
    "    fill_color=None, size = 9, color=\"darkorange\",\n",
    "    tools=[detHoverTool])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we overlay the detected sources on the image. The `*` operator  is used to overlay one Element on to another.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset the tools on the image and add a hover on the detections.\n",
    "rasterize(img).opts(tools=[]) * detections.opts(tools=[detHoverTool])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now  mouse-over the sources and get the coordinates of the detections. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.0 Interactive Image Exploration with with Holoviews Streams and DynamicMap\n",
    "\n",
    "Now let's add some interactive exploration capability using Holoviews [Streams](http://holoviews.org/user_guide/Streaming_Data.html) and [DynamicMap](https://holoviews.org/reference/containers/bokeh/DynamicMap.html). A DynamicMap is an explorable multi-dimensional wrapper around a callable that returns HoloViews objects. The core concept behind a stream is simple: it defines one or more parameters that can change over time that automatically refreshes code depending on those parameter values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First create a DynamicMap with a box stream so that we can explore selected sections of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundsxy = (0, 0, 0, 0)\n",
    "box = streams.BoundsXY(source=img, bounds=boundsxy)\n",
    "dynamicMap = hv.DynamicMap(lambda bounds: hv.Bounds(bounds), streams=[box])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the image and overlay the DynamicMap\n",
    "rasterize(img) * dynamicMap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the interactive callback features on the image plots, such as the selection box (the icon of the box with a '+' in the lower right corner), we can explore regions of the image.  Use the box select tool on the image above to select a region and then execute the cell below to get the box boundary coordinates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is another version of the image with a [tap stream](http://holoviews.org/reference/streams/bokeh/Tap.html) instead of box select. A Tap stream allows you to click or 'tap' a position to interact with a plot. Try zooming in on an interesting part of the image and then 'tap' somewhere to place an 'X' marker. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posxy = hv.streams.Tap(source=img, x=0.5 * calexp.getDimensions()[0],\n",
    "                       y=0.5 * calexp.getDimensions()[1])\n",
    "marker = hv.DynamicMap(lambda x, y: hv.Points([(x, y)]), streams=[posxy])\n",
    "rasterize(img)* marker.opts(color='white', marker='x', size=20)"
   ]
  },
  {
   "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 scaled/raw value at position (%.3f, %.3f) is %.3f / %.3f' %\n",
    "      (posxy.x, posxy.y, scaledImage[-int(posxy.y), int(posxy.x)], \n",
    "       calexp.image.array[-int(posxy.y), int(posxy.x)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.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",
    " 2. In the image display sections, try using the coadd image instead of the calexp image. \n",
    " \n",
    " 3. In section 2.3, try extracting additional information about the detected sources from the calexpSrc table and adding it to the custom hover tool. For example, the corresponding RA/DEC or the PSF flux. \n",
    " \n",
    " 3. Try using a different stream function to interact with the images in section 3  "
   ]
  }
 ],
 "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
}