03_Image_Display_and_Manipulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<b>Displaying images using the LSST astronomical framework (afw) library</b> <br>\n",
    "Last verified to run on 2021-10-12 with LSST Science Pipelines release w_2021_40 <br>\n",
    "Contact authors: Alex Drlica-Wagner, Jeff Carlin <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "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 is based on the [`AFW_Display_Demo.ipynb`](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Visualization/AFW_Display_Demo.ipynb) notebook originally written by Brant Robertson and maintained by the LSST Stack Club. Please consider acknowledging Brant Robertson, Alex Drlica-Wagner, and Jeff Carlin in any publications or software releases that make use of this notebook's contents.\n",
    "\n",
    "More examples of the use of `lsst.afw.display` can be found in the [LSST Science Pipelines documentation](https://pipelines.lsst.io/getting-started/display.html). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\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.butler` 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",
    "### Set Up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the Stack are we using?\n",
    "! echo $IMAGE_DESCRIPTION\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### 1. 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 below, including the `matplotlib.pyplot` plotting sublibrary and the [`astropy.wcs`](https://docs.astropy.org/en/stable/wcs/index.html) package for dealing with World Coordinate Systems (WCS). 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 matplotlib.pyplot as plt      # imports matplotlib.pyplot as plt\n",
    "import warnings                      # imports the warnings library\n",
    "import gc                            # imports python's garbage collector\n",
    "from astropy.wcs import WCS          # imports astropy's World Coordinate System function WCS"
   ]
  },
  {
   "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": [
    "# prevent some helpful but ancillary warning messages from printing\n",
    "#   during some LSST DM Release calls\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matplotlib stores the data array associated with an image that is plotted. Since the LSST images are large (~4k x 4k pixels), this can eventually lead to a memory overflow, which will cause the notebook kernel to die. To mitigate this issue, we define a function to clean up after we plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_figure(fig):\n",
    "    \"\"\"Remove a figure to reduce memory footprint. \"\"\"\n",
    "    # get the axes and clear their images\n",
    "    for ax in fig.get_axes():\n",
    "        for im in ax.get_images():\n",
    "            im.remove()\n",
    "    fig.clf()      # clear the figure\n",
    "    plt.close(fig) # close the figure\n",
    "    gc.collect()   # call the garbage collector"
   ]
  },
  {
   "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, and set some other plotting parameters to make things look nice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up some plotting defaults:\n",
    "\n",
    "params = {'axes.labelsize': 28,\n",
    "          'font.size': 24,\n",
    "          'legend.fontsize': 14,\n",
    "          'xtick.major.width': 3,\n",
    "          'xtick.minor.width': 2,\n",
    "          'xtick.major.size': 12,\n",
    "          'xtick.minor.size': 6,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 3,\n",
    "          'axes.linewidth': 3,\n",
    "          'axes.labelweight': 3,\n",
    "          'axes.titleweight': 3,\n",
    "          'ytick.major.width': 3,\n",
    "          'ytick.minor.width': 2,\n",
    "          'ytick.major.size': 12,\n",
    "          'ytick.minor.size': 6,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [8, 8],\n",
    "          'figure.facecolor': 'White'\n",
    "          }\n",
    "\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Load the LSST Science Pipelines\n",
    "\n",
    "First, we load the `lsst.afw.display` library to gain access to the image visualization routines we'd like to use, and the `lsst.daf.butler` library, which is used to access data products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load lsst.afw.display to gain access to image visualization routines.\n",
    "import lsst.afw.display as afwDisplay\n",
    "# load the Butler, which provides programmatic access to LSST data products.\n",
    "from lsst.daf.butler import Butler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Load the Data to Visualize\n",
    "\n",
    "To display an image, we must first load some data. These data have been processed with the LSST Science Pipelines, and are organized in a structure that enables us to access them through the `Butler`. For more information on the `Butler`, see [lsst.daf.butler](https://pipelines.lsst.io/modules/lsst.daf.butler/index.html).\n",
    "\n",
    "The DP0.1 data set contains simulated images from the LSST DESC Data Challenge 2 (DC2). These data are available in an S3 bucket `s3://butler-us-central1-dp01`. We access a single image from a specific visit (`192350`) and detector (`175`).\n",
    "\n",
    "Once we define a string that contains the data directory, we start the `Butler` instance using the `lsst.daf.butler` library 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 (\"band\"), CCD (\"detector\", and visit we wish to view. To do this, we define a dictionary with the required information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = {'visit': 192350, 'detector': 175, 'band': 'i'}\n",
    "# Note: because the combination of visit+detector already uniquely identifies\n",
    "# the exposure, specifying \"band\" above is unnecessary.\n",
    "\n",
    "repo = 's3://butler-us-central1-dp01'\n",
    "collection = '2.2i/runs/DP0.1'\n",
    "butler = 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": [
    "### 4.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 use `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": [
    "# Use lsst.afw.display with the matplotlib backend\n",
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "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 `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. To 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": [
    "# create a matplotlib.pyplot figure\n",
    "fig, ax = plt.subplots()\n",
    "# get an alias to the lsst.afw.display.Display() method\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "# set the image stretch algorithm and range\n",
    "display.scale('asinh', 'zscale')\n",
    "# load the image into the display\n",
    "display.mtv(calexp.image)\n",
    "# show the corresponding pyplot figure\n",
    "plt.show()\n",
    "# clean up memory\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the image axes in sky coordinates instead of pixel coordinates, a simple option is to use astropy's World Coordinate System (WCS) package, along with matplotlib.pyplot's ``subplot``, ``imshow``, and ``grid`` functions. \n",
    "Recall that we imported ``matplotlib.pyplot`` as ``plt`` already, and that we imported the ``astropy.wcs.WCS`` function as simply ``WCS``.\n",
    "Find more information about [imshow](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.imshow.html) and [colormaps](https://matplotlib.org/stable/tutorials/colors/colormaps.html) (``cmap``)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "# Set the figure's projection to be the WCS of the calexp\n",
    "plt.subplot(projection=WCS(calexp.getWcs().getFitsMetadata()))\n",
    "# Display the calexp image data array using the gray colormap (cmap)\n",
    "#  and use approximately the same min and max scale values as above\n",
    "im = plt.imshow(calexp.image.array, cmap='gray', vmin=-200.0, vmax=400, origin='lower')\n",
    "# Add solid white grid lines\n",
    "plt.grid(color='white', ls='solid')\n",
    "# Label axes\n",
    "plt.xlabel('Right Ascension')\n",
    "plt.ylabel('Declination')\n",
    "plt.show()\n",
    "# Clean up memory\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often you may want to plot two images side-by-side. This can be done by creating matplotlib subplots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14, 7))\n",
    "plt.sca(ax[0])  # set the first axis as current\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "display1.scale('linear', 'zscale')\n",
    "display1.mtv(calexp.image)\n",
    "plt.sca(ax[1])  # set the second axis as current\n",
    "display2 = afwDisplay.Display(frame=fig)\n",
    "display2.mtv(calexp.mask)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to plot the mask on top of the image using the `calexp.maskedImage`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('linear', 'zscale')\n",
    "display.mtv(calexp.maskedImage)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Congratulations!** We've plotted an image using `lsst.afw.display`!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.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 Stack Club [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\n",
    "   (as a percentage: 0 = opaque, 100 = transparent)\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": [
    "# create a matplotlib.pyplot figure\n",
    "fig, ax = plt.subplots()\n",
    "# get an alias to the lsst.afw.display.Display() method\n",
    "afw_display = afwDisplay.Display(frame=fig)\n",
    "# set the image stretch algorithm and range\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "# set the transparency of the mask plane (0%=opaque, 100%=transparent)\n",
    "afw_display.setMaskTransparency(80)\n",
    "# set the color for a single plane in the mask\n",
    "afw_display.setMaskPlaneColor('DETECTED', 'blue')\n",
    "# load the image and mask plane into the display\n",
    "afw_display.mtv(calexp)\n",
    "# show the corresponding pyplot figure\n",
    "plt.show()\n",
    "# clean up memory\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `afw_display` object contains more information about the mask planes that can be accessed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the colors associated to each plane in the mask\n",
    "print(\"Mask plane bit definitions:\\n\", afw_display.getMaskPlaneColor())\n",
    "print(\"\\nMask plane methods:\\n\")\n",
    "help(afw_display.setMaskPlaneColor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. 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 = [fun for fun in dir(display) if callable(getattr(display, fun))]\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": [
    "### 6. Extract a cutout image\n",
    "\n",
    "Say we want to grab a cutout of the DP0.1 coadded images at a specific location. In order to do this, we need a few other packages from the LSST Science Pipelines. In particular, we'll need access to the geometry and coordinate packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.geom as geom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will define a convenience function to extract a cutout from a deep coadd image, if given an RA, Dec position and desired image size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutout_coadd(butler, ra, dec, band='r', datasetType='deepCoadd',\n",
    "                 skymap=None, cutoutSideLength=51, **kwargs):\n",
    "    \"\"\"\n",
    "    Produce a cutout from a coadd at the given ra, dec position.\n",
    "\n",
    "    Adapted from DC2 tutorial notebook by Michael Wood-Vasey.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Servant providing access to a data repository\n",
    "    ra: float\n",
    "        Right ascension of the center of the cutout, in degrees\n",
    "    dec: float\n",
    "        Declination of the center of the cutout, in degrees\n",
    "    band: string\n",
    "        Filter of the image to load\n",
    "    datasetType: string ['deepCoadd']\n",
    "        Which type of coadd to load.  Doesn't support 'calexp'\n",
    "    skymap: lsst.afw.skyMap.SkyMap [optional]\n",
    "        Pass in to avoid the Butler read.  Useful if you have lots of them.\n",
    "    cutoutSideLength: float [optional]\n",
    "        Size of the cutout region in pixels.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \"\"\"\n",
    "    radec = geom.SpherePoint(ra, dec, geom.degrees)\n",
    "    cutoutSize = geom.ExtentI(cutoutSideLength, cutoutSideLength)\n",
    "\n",
    "    if skymap is None:\n",
    "        skymap = butler.get(\"skyMap\")\n",
    "\n",
    "    # Look up the tract, patch for the RA, Dec\n",
    "    tractInfo = skymap.findTract(radec)\n",
    "    patchInfo = tractInfo.findPatch(radec)\n",
    "    xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))\n",
    "    bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)\n",
    "    patch = tractInfo.getSequentialPatchIndex(patchInfo)\n",
    "\n",
    "    coaddId = {'tract': tractInfo.getId(), 'patch': patch, 'band': band}\n",
    "    parameters = {'bbox': bbox}\n",
    "\n",
    "    cutout_image = butler.get(datasetType, parameters=parameters,\n",
    "                              immediate=True, dataId=coaddId)\n",
    "\n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we extract a cutout image, let's take a look at what a full 4k x 4k pixel \"patch\" image looks like. Previously we were extracting a single visit image spanning only one detector. For the rest of this notebook, we will look at coadded images made up of multiple exposures that have been combined.\n",
    "\n",
    "Let's grab a \"deepCoadd_calexp\" image and display the entire patch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = {'tract': 4431, 'patch': 17, 'band': 'i'}\n",
    "\n",
    "# Retrieve the data using the `butler` instance and its function `get()`\n",
    "coadd_calexp = butler.get('deepCoadd_calexp', **dataId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a matplotlib.pyplot figure\n",
    "fig, ax = plt.subplots()\n",
    "# get an alias to the lsst.afw.display.Display() method\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "# set the image stretch algorithm and range\n",
    "display.scale('asinh', 'zscale')\n",
    "# load the image into the display\n",
    "display.mtv(coadd_calexp.image)\n",
    "plt.show()\n",
    "# clean up memory\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow - check out that rich galaxy cluster in the lower-left portion of the image! Let's use the function we defined above to get a cutout around the position of that cluster.\n",
    "\n",
    "(You may have noticed, though, that the image above is displaying pixel numbers rather than RA, Dec coordinates. As a quick aside, we'll show you how to get the sky coordinates corresponding to an XY pixel position.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we need to extract the WCS solution, which provides the mapping\n",
    "#   between XY pixel values and sky coordinates:\n",
    "wcs = coadd_calexp.getWcs()\n",
    "\n",
    "# Print the WCS info to see what it contains:\n",
    "print(wcs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The cluster seems to be centered at about (X, Y) = (12500, 8500).\n",
    "# We can use the \"pixelToSky\" method of the WCS to get the sky coordinates:\n",
    "radec = wcs.pixelToSky(12500, 8500)\n",
    "\n",
    "print(radec.getRa().asDegrees(), radec.getDec().asDegrees())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use this RA, Dec position to extract a cutout image that is centered on the galaxy cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a position at roughly the center of the galaxy cluster:\n",
    "ra, dec = radec.getRa().asDegrees(), radec.getDec().asDegrees()\n",
    "cutout_image = cutout_coadd(butler, ra, dec, datasetType='deepCoadd',\n",
    "                            cutoutSideLength=501)\n",
    "print(\"The size of the cutout in pixels is: \", cutout_image.image.array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "afw_display = afwDisplay.Display(frame=fig)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(cutout_image.image)\n",
    "ax.axis('off')\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Plot an RGB composite image\n",
    "\n",
    "That's a pretty cluster, but what if we really want to know about the colors of the stars and galaxies in that image? To do that, we can extract images taken with three different filters, then assign those images to the RGB channels of a color image.\n",
    "\n",
    "First we'll load a tool from Astropy that helps create RGB images (see [this Astropy guide to creating RGB images](https://docs.astropy.org/en/stable/visualization/rgb.html?highlight=make_lupton_rgb) for more info), and make a convenience function to combine the steps of making an RGB image. (This function was originally written by Aaron Watkins and Nushkia Chamba for their [Stack Club Course](https://github.com/LSSTScienceCollaborations/StackClubCourse) project exploring surface photometry in the LSST Science Pipelines.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.visualization import make_lupton_rgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createRGB(image, bgr=\"gri\", stretch=1, Q=10, scale=None):\n",
    "    \"\"\"\n",
    "    Create an RGB color composite image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    image : `MultibandExposure`\n",
    "        `MultibandExposure` to display.\n",
    "    bgr : sequence\n",
    "        A 3-element sequence of filter names (i.e., keys of the exps dict)\n",
    "        indicating what band to use for each channel. If `image` only has\n",
    "        three filters then this parameter is ignored and the filters\n",
    "        in the image are used.\n",
    "    stretch: int\n",
    "        The linear stretch of the image.\n",
    "    Q: int\n",
    "        The Asinh softening parameter.\n",
    "    scale: list of 3 floats, each less than 1. (default: None)\n",
    "        Re-scales the RGB channels.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    rgb: ndarray\n",
    "        RGB (integer, 8-bits per channel) colour image as an NxNx3 numpy array.\n",
    "    \"\"\"\n",
    "\n",
    "    # If the image only has 3 bands, reverse the order of the bands\n",
    "    #   to produce the RGB image\n",
    "    if len(image) == 3:\n",
    "        bgr = image.filters\n",
    "\n",
    "    # Extract the primary image component of each Exposure with the\n",
    "    #   .image property, and use .array to get a NumPy array view.\n",
    "\n",
    "    if scale is None:\n",
    "        r_im = image[bgr[2]].array  # numpy array for the r channel\n",
    "        g_im = image[bgr[1]].array  # numpy array for the g channel\n",
    "        b_im = image[bgr[0]].array  # numpy array for the b channel\n",
    "    else:\n",
    "        # manually re-scaling the images here\n",
    "        r_im = image[bgr[2]].array * scale[0]\n",
    "        g_im = image[bgr[1]].array * scale[1]\n",
    "        b_im = image[bgr[0]].array * scale[2]\n",
    "\n",
    "    rgb = make_lupton_rgb(image_r=r_im,\n",
    "                          image_g=g_im,\n",
    "                          image_b=b_im,\n",
    "                          stretch=stretch, Q=Q)\n",
    "    # \"stretch\" and \"Q\" are parameters to stretch and scale the pixel values\n",
    "\n",
    "    return rgb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the `MultiBandExposure` object type from the LSST Science Pipelines `afw.image` package. One could also create a similar function that simply accepts three separate images without combining them into a MultiBand object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.image import MultibandExposure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutout_image_g = cutout_coadd(butler, ra, dec, band='g',\n",
    "                              datasetType='deepCoadd', cutoutSideLength=701)\n",
    "cutout_image_r = cutout_coadd(butler, ra, dec, band='r',\n",
    "                              datasetType='deepCoadd', cutoutSideLength=701)\n",
    "cutout_image_i = cutout_coadd(butler, ra, dec, band='i',\n",
    "                              datasetType='deepCoadd', cutoutSideLength=701)\n",
    "\n",
    "# Multiband exposures need a list of images and filters\n",
    "coadds = [cutout_image_g, cutout_image_r, cutout_image_i]\n",
    "coadds = MultibandExposure.fromExposures(['g', 'r', 'i'], coadds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to create the RGB (where, in case you're unfamiliar, the RGB denotes \"red,\" \"green,\" and \"blue\" channels) image. We will assign the g, r, and i-bands to the B, G, and R channels (i.e., g-band, as the bluest filter, will be in the blue channel, r-band in the green channel, and i-band as red). \n",
    "\n",
    "Finally, we will plot two versions of the image to demonstrate how you can change the relative scaling of each RGB channel to produce an image better highlights certain features. You can experiment with the \"scale = []\" keyword, or including different bands besides _gri_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(20, 20), nrows=1, ncols=2)\n",
    "\n",
    "# original make_lupton_rgb without any scaling\n",
    "rgb_original = createRGB(coadds.image, bgr=['g', 'r', 'i'], scale=None)\n",
    "ax[0].imshow(rgb_original, origin='lower')\n",
    "ax[0].set_title('original', fontsize=30)\n",
    "\n",
    "# make_lupton_rgb with scaled rgb channels\n",
    "ax[1].set_title('re-scaled', fontsize=30)\n",
    "rgb_scaled = createRGB(coadds.image, bgr=['g', 'r', 'i'],\n",
    "                       scale=[0.6, 0.7, 1.0])\n",
    "ax[1].imshow(rgb_scaled, origin='lower')\n",
    "\n",
    "ax[0].set_axis_off()\n",
    "ax[1].set_axis_off()\n",
    "plt.show()\n",
    "\n",
    "# clean up memory\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations - you have now learned how to display images, create cutouts centered on a given position, and make multi-band color images. Enjoy exploring the DP0 images!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additional 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
}