Visualization/AFW_Display_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# **Displaying images using astronomical framework library (lsst.afw.display)**\n",
    "\n",
    "<br>Author(s): **Brant Robertson** ([@brantr](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@brantr))\n",
    "<br>Maintainer(s): **Alex Drlica-Wagner** ([@kadrlica](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica))\n",
    "<br>Level: **Introductory**\n",
    "<br>Last Verified to Run: **2021-09-11**\n",
    "<br>Verified Stack Release: **w_2021_33**\n",
    "\n",
    "### Learning Objectives:\n",
    "\n",
    "In this tutorial we will \n",
    "\n",
    "1. Show how to access the `lsst.afw.display` routines.\n",
    "2. Use the LSST data Butler to access processed data and inspect it visually.\n",
    "\n",
    "This tutorial is designed to help users get a brief feel for the `lsst.afw.display` library that enables the visual inspection of data. The [`lsst.afw` library](https://github.com/lsst/afw) provides an \"Astronomical Framework\" (afw) while the `lsst.daf.*` libraries (see, e.g., [daf_base](https://github.com/lsst/daf_base)) provides a Data Access Framework (daf). Both libraries are used in this tutorial, with the `lsst.daf.persistence` library used to access a calibrated exposure (calexp) and the `lsst.afw.display` library used to show the exposure image on the screen.\n",
    "\n",
    "This tutorial made use of the [`LowSurfaceBrightness.ipynb` StackClub notebook](https://nbviewer.jupyter.org/github/LSSTScienceCollaborations/StackClub/blob/rendered/SourceDetection/LowSurfaceBrightness.nbconvert.ipynb) by [Alex Drlica-Wagner](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica). More examples of the use of `lsst.afw.display` can be found in the [Stack ](https://pipelines.lsst.io/getting-started/display.html).\n",
    "\n",
    "### Logistics\n",
    "This notebook is intended to be run on a <span style=\"color:blue;font-weight:bold\">LARGE</span> instance of the RSP at `lsst-lsp-stable.ncsa.illinois.edu` or `data.lsst.cloud` from a local git clone of the [StackClub](https://github.com/LSSTScienceCollaborations/StackClub) repo.\n",
    "\n",
    "\n",
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Site, host, and stack version\n",
    "! echo $EXTERNAL_INSTANCE_URL\n",
    "! echo $HOSTNAME\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Step 0) Import Common Python Libraries\n",
    "\n",
    "The [`matplotlib`](https://matplotlib.org/), [`numpy`](http://www.numpy.org/), and [`astropy`](http://www.astropy.org/) libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis. We will use these packages in common ways below, including the `matplotlib.pyplot` plotting sublibrary. We also import the [`warnings` library](https://docs.python.org/2/library/warnings.html) to prevent some routine warning messages from printing to the screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#allow for matplotlib to create inline plots in our notebook\n",
    "%matplotlib inline                   \n",
    "import os                            #imports os library\n",
    "import numpy as np                   #imports numpy with the alias np\n",
    "import matplotlib.pyplot as plt      #imports matplotlib.pyplot as plt\n",
    "import warnings                      #imports the warnings library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go ahead and import from `astropy` the image stretch limits from the familiar `zscale()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.visualization import ZScaleInterval     #This function allows use to use the `zscale()` rescaling limits function familiar from, e.g., DS9, to adjust the image stretch.\n",
    "zscale = ZScaleInterval()                            #create an alias to the `ZScaleInterval()` function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let the kernel know that we're happy not to have some useful warnings printed during this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "warnings.simplefilter(\"ignore\", category=FutureWarning)   #prevent some helpful but ancillary warning messages from printing during some LSST DM Release calls\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)     #prevent some helpful but ancillary warning messages from printing during some LSST DM Release calls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a last preparatory task, we set the parameters of `matplotlib.pyplot` to give us a large default size for an image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8.0, 8.0)               #set a large default size for our images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1) Loading the LSST DM Stack\n",
    "\n",
    "First, we load some LSST stack libraries to gain access to the image data and visualization routines we'd like to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.daf.butler      as dafButler      #load the Butler to access data \n",
    "import lsst.afw.display     as afwDisplay     #load lsst.afw.display to gain access to image visualization routines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2) Importing Data to Visualize\n",
    "\n",
    "To display an image, we must first load some data. These data have been processed with the LSST DM Stack, and are organized in a structure that enables us to access them through the DM Stack `Butler`. For more information on the `Butler`, see [lsst.daf.butler](https://pipelines.lsst.io/modules/lsst.daf.butler/index.html).\n",
    "\n",
    "In this tutorial, we provide access to two different data sets:\n",
    "1. The first data set contains real imaging data from a reprocessing of `HSC` data available at NCSA in the data directory `/datasets/hsc/repo/rerun/RC/v20_0_0_rc1/DM-25349`. We access an image from a specific visit (`38938`) and ccd (`32`). This happens to be an HSC z-band exposure.\n",
    "2. The second data set contains simulated images from the LSST DESC Data Challenge 2 (DC2). These data are available at NCSA in the data directory `/datasets/DC2/DR6/Run2.2i/patched/2021-02-10/rerun/run2.2i-calexp-v1/`. We access a single image from a specific visit (`512055`), raft (`R20`), and detector (`76`). This happens to be an i-band exposure.\n",
    "\n",
    "Once we define a string that contains the data directory, we start the `Butler` instance using the `lsst.daf.persistence` library alias `dafPersist` and its `Butler` class. The `Butler` object is initialized with a string containing the data directory we wish to access. Running the cell may take a few moments.\n",
    "\n",
    "With the `Butler` instance now generated using our data directory, we can retrieve the desired calibrated exposure by telling the butler which filter, CCD, and visit we wish to view. To do this, we definie dictionary with the required information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Location of the DC2 data repository on this site\n",
    "URL = os.getenv('EXTERNAL_INSTANCE_URL')\n",
    "if URL.endswith('data.lsst.cloud'): # IDF\n",
    "    repo = \"s3://butler-us-central1-dp01\"\n",
    "elif URL.endswith('ncsa.illinois.edu'): # NCSA\n",
    "    repo = \"/repo/dc2\"\n",
    "else:\n",
    "    raise Exception(f\"Unrecognized URL: {URL}\")\n",
    "\n",
    "dataset='DC2'\n",
    "collection='2.2i/runs/DP0.1'\n",
    "dataId = {'visit': 227982, 'raftName': 'R31', 'detector': 129}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the butler\n",
    "butler = dafButler.Butler(repo,collections=collection)\n",
    "\n",
    "# Retrieve the data using the `butler` instance and its function `get()`\n",
    "calexp = butler.get('calexp', **dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3.1) Use AFWDisplay to Visualize the Image\n",
    "\n",
    "Now, with a `Butler` instance defined and a calibrated exposure retrieved, we can use [`lsst.afw.display`](https://github.com/lsst/afw/tree/master/python/lsst/afw/display) to visualize the data.  The next task is to let AFWDisplay know that we want it to enroll `matplotlib` as our default display backend. To do this, we use the `setDefaultBackend()` function. Remember that we made an alias to `lsst.afw.display` called `afwDisplay`, so we'll use that to call `setDefaultBackend()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')  # Use lsst.afw.display with the matplotlib backend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now set to display the image. To do this, we:\n",
    "\n",
    "1. First create a `matplotlib.pyplot` figure using `plt.figure()` -- this will be familiar to anyone with experience using `matplotlib`.\n",
    "2. Then create an alias to the `lsst.afw.display.Display` method that will allow us to display the data to the screen.  This alias will be called `afw_display`.\n",
    "3. Before showing the data on the screen, we have to decide how to apply an image stretch given the data. The algorithm we'll use is `asinh` familiar from SDSS images, with a range of values set by `zscale`. To do this, we use the `scale()` function provided by `lsst.afw.display`. See the `scale()` function definition in the [`interface.py` file of the lsst.afw.display library](https://github.com/lsst/afw/blob/master/python/lsst/afw/display/interface.py).\n",
    "4. Finally, we can display the image. Do do this, we provide the `mtv()` method the `image` member of our calibrated image retrieved by the `butler`. We can then use `plt.show()` to display our figure.\n",
    "\n",
    "All these tasks are best done within the same notebook cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()                        #create a matplotlib.pyplot figure\n",
    "display = afwDisplay.Display()      #get an alias to the lsst.afw.display.Display() method\n",
    "display.scale('asinh', 'zscale')    #set the image stretch algorithm and range\n",
    "display.mtv(calexp.image)           #load the image into the display\n",
    "plt.show()                          #show the corresponding pyplot figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often you may want to plot two images side-by-side. This can be done by creating matplotlib subplots. Here we plot the image and its mask side by side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,2,figsize=(14,7))\n",
    "\n",
    "plt.sca(ax[0]) # set the first axes as current\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "display1.scale('linear', 'zscale')\n",
    "display1.mtv(calexp.image)\n",
    "\n",
    "plt.sca(ax[1]) # set the second axes as current\n",
    "display2 = afwDisplay.Display(frame=fig)\n",
    "display2.mtv(calexp.mask)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to plot the mask on top of the image using the `calexp.maskedImage`. In this case, the image pixel values are plotted in grayscale with the mask values are overplotted in semi-transparent color. We explort the mask plane a bit more in the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('linear', 'zscale')\n",
    "display.mtv(calexp.maskedImage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Congrats!** We've plotted an image using `lsst.afw.display`!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3.2) Use AFWDisplay to Visualize the Image and Mask Plane\n",
    "\n",
    "The `calexp` returned by the butler contains more than just the image pixel values (see the [calexp tutorial](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Basics/Calexp_guided_tour.ipynb) for more details). One other component is the mask plane associated with the image. `AFWDisplay` provides a nice pre-packaged interface for overplotting the mask associated with an image. A mask is composed of a set of \"mask planes\", 2D binary bit maps corresponding to pixels that are masked for various reasons (see [here](https://pipelines.lsst.io/v/DM-11392/getting-started/display.html#interpreting-displayed-mask-colors) for more details)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll follow the same steps as above to display the image, but we'll add a few modifications\n",
    "\n",
    "1. We explicitly set the transparency of the overplotted mask (0 = transparent, 1 = opaque)\n",
    "2. We explicitly set the color of the 'DETECTED' mask plane to 'blue' (i.e. all pixels associated with detected objects).\n",
    "3. We pass the full `calexp` object to `mtv` instead of just the image plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()                                     #create a matplotlib.pyplot figure\n",
    "display = afwDisplay.Display()               #get an alias to the lsst.afw.display.Display() method\n",
    "display.scale('asinh', 'zscale')             #set the image stretch algorithm and range\n",
    "display.setMaskTransparency(0)               #set the transparency of the mask plane (0=opaque; 100=transparent)\n",
    "display.setMaskPlaneColor('DETECTED','pink') #set the color for a single plane in the mask\n",
    "display.mtv(calexp)                          #load the image and mask plane into the display\n",
    "plt.show()                                       #show the corresponding pyplot figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `display` object contains more information about the mask planes that can be accessed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Mask plane bit definitions:\\n\", display.getMaskPlaneColor())  # Print the colors associated to each plane in the mask\n",
    "print(\"\\nMask plane methods:\\n\")\n",
    "help(display.setMaskPlaneColor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A compact way to print some mask information\n",
    "for maskName, maskBit in calexp.mask.getMaskPlaneDict().items(): \n",
    "    print(f'{maskName:18s} (bit=2**{maskBit:02d}): {display.getMaskPlaneColor(maskName)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that some mask plane bits have their colors unset (i.e., `None`). These mask planes have colors assigned on the fly at draw-time by `mtv`. Unfortunately, there isn't a clear way to find out what color was assigned, and this color may change if additional mask planes are added or removed. The best way to have reproducible mask plane colors is to assign an explicit color to each mask plane with `setMaskPlaneColor`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3.3) Use AFWDisplay to Visualize a Coadd Image\n",
    "\n",
    "Now we are going to repeat some of the plotting using a coadd data product. First we grab the coadd image and then display it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First use the Butler to get the coadd\n",
    "coaddId = {'band':'i', 'skymap': 'DC2', 'tract': 4851, 'patch': 29}\n",
    "coadd = butler.get('deepCoadd', **coaddId)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the coadd; look at all the mask bits!\n",
    "plt.figure()                          \n",
    "display = afwDisplay.Display()    \n",
    "display.scale('asinh', 'zscale')  \n",
    "display.mtv(coadd)                \n",
    "plt.show()                            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow, look at that mask plane! We can investigate the mask bits as we did for the `calexp`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A compact way to print some mask information\n",
    "print(\"Mask Plane Dict\")\n",
    "for maskName, maskBit in coadd.mask.getMaskPlaneDict().items(): \n",
    "    print(f'{maskName:18s} (bit=2**{maskBit:02d}): {display.getMaskPlaneColor(maskName)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that much of the mask plane is coming from the edges of the sensors that go into building the coadd. To remove this mask bit, we make the `SENSOR_EDGE` and `INEXACT_PSF` mask planes transparent. The remaining streaky features in the mask plane coming from the bright objects are coming from the `REJECTED` bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's make the SENSOR_EDGE mask plane transparent\n",
    "plt.figure()                                     \n",
    "display = afwDisplay.Display()               \n",
    "display.scale('asinh', 'zscale') \n",
    "display.setMaskTransparency(100,name='SENSOR_EDGE')\n",
    "display.setMaskTransparency(100,name='INEXACT_PSF')\n",
    "#display.setMaskTransparency(100,name='REJECTED')\n",
    "display.mtv(coadd)\n",
    "plt.show()                                       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also pan and zoom\n",
    "plt.figure()                                     \n",
    "display = afwDisplay.Display()               \n",
    "display.scale('asinh', 'zscale') \n",
    "display.mtv(coadd.image)\n",
    "display.zoom(16)\n",
    "display.pan(5150, 18650)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4) More Information about lsst.afw.display\n",
    "\n",
    "To get some more information about `lsst.afw.display`, we can print the method list to see what's available. The next cell will print `lsst.afw.display` methods to the screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "method_list = [func for func in dir(display) if callable(getattr(display, func))]\n",
    "print(method_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you'd like to learn more about any given function, please see the [`lsst.afw.display` source code](https://github.com/lsst/afw/tree/master/python/lsst/afw/display).\n",
    "\n",
    "You can also read the API documentation about the above functions using the Jupyter notebook `help()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(display.scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(display.mtv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Further Documentation**\n",
    "\n",
    "If you'd like some more information on `lsst.afw.display`, please have a look at the following websites:\n",
    "\n",
    "* [Info on image indexing conventions.](https://pipelines.lsst.io/modules/lsst.afw.image/indexing-conventions.html)  \n",
    "* [afw.display Doxygen website](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/namespacelsst_1_1afw_1_1display.html)  \n",
    "* [afw.display GitHub website](https://github.com/RobertLuptonTheGood/afw/tree/master/python/lsst/afw/display)  \n",
    "* [Getting Started on Image Display (pipelines.lsst.io)](https://pipelines.lsst.io/getting-started/display.html)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.8.8"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}