DP02_03a_Image_Display_and_Manipulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Display and Manipulation\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br>\n",
    "Contact authors: Alex Drlica-Wagner, Jeff Carlin <br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_50 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Description:** Learn to display and manipulate images using the LSST Science Pipelines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Skills:** Display and manipulate images, explore image mask planes, create cutouts and RGB images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**LSST Data Products:** Butler calexp and deepCoadd images, coadd objectTable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Packages:** lsst.afw.display, lsst.daf.butler, lsst.geom, lsst.afw.image"
   ]
  },
  {
   "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.\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": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>.\n",
    "Questions are welcome as new topics in the <a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> of the Rubin Community Forum. Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction\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 image 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)) provide a \"Data Access Framework\" (daf). Both libraries are used in this tutorial, with the `lsst.daf.butler` library used to access image data and the `lsst.afw.display` library used to show the exposure image on the screen using the matplotlib backend."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### 1.1. Package imports\n",
    "\n",
    "Below, the `matplotlib.pyplot` sublibrary is imported for plotting, and the [`astropy.wcs`](https://docs.astropy.org/en/stable/wcs/index.html) package is imported for dealing with World Coordinate Systems (WCS).\n",
    "The [`matplotlib`](https://matplotlib.org/) and [`astropy`](http://www.astropy.org/) libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis.\n",
    "The `gc` (garbage collector) package is imported to help clear the memory of large plots.\n",
    "\n",
    "The `lsst.afw.display` and `lsst.afw.image` modules are loaded to gain access to the image visualization routines we'd like to use, and the `lsst.daf.butler` and `lsst.rsp.get_tap_service` modules are used to access data products.\n",
    "The `lsst.geom` package is used for dealing with sky coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from astropy.wcs import WCS\n",
    "from astropy.visualization import make_lupton_rgb\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "import gc\n",
    "\n",
    "import lsst.afw.display as afwDisplay\n",
    "from lsst.afw.image import MultibandExposure\n",
    "from lsst.daf.butler import Butler\n",
    "from lsst.rsp import get_tap_service\n",
    "import lsst.geom as geom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2.1. Functions\n",
    "\n",
    "Matplotlib stores the data array associated with an image that is plotted. Since the LSST Charge-Coupled Device (CCD) detector 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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def remove_figure(fig):\n",
    "    \"\"\"\n",
    "    Remove a figure to reduce memory footprint.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fig: matplotlib.figure.Figure\n",
    "        Figure to be removed.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    None\n",
    "    \"\"\"\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": [
    "<a id=\"cutout_coadd\"></a>\n",
    "\n",
    "Next we define a function to produce a cutout image from a coadd at the user-provided RA, Dec position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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",
    "        Helper object 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",
    "                              dataId=coaddId)\n",
    "\n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function makes cutouts from `calexp` images instead of `deepCoadd` images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def cutout_calexp(butler, ra, dec, visit, detector, cutoutSideLength=51, **kwargs):\n",
    "    \n",
    "    \"\"\"\n",
    "    Produce a cutout from a calexp at the given ra, dec position.\n",
    "\n",
    "    Adapted from cutout_coadd which was adapted from a DC2 tutorial\n",
    "    notebook by Michael Wood-Vasey.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Helper object 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",
    "    visit: int\n",
    "        Visit id of the calexp's visit\n",
    "    detector: int\n",
    "        Detector for the calexp\n",
    "    cutoutSideLength: float [optional]\n",
    "        Size of the cutout region in pixels.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \"\"\"\n",
    "    \n",
    "    dataId = {'visit': visit, 'detector': detector}    \n",
    "    radec = geom.SpherePoint(ra, dec, geom.degrees)\n",
    "    cutoutSize = geom.ExtentI(cutoutSideLength, cutoutSideLength)    \n",
    "    calexp_wcs = butler.get('calexp.wcs', **dataId)\n",
    "    xy = geom.PointI(calexp_wcs.skyToPixel(radec))\n",
    "    bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)\n",
    "    parameters = {'bbox': bbox}\n",
    "    cutout_image = butler.get('calexp', parameters=parameters, **dataId)\n",
    "\n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"create_rgb\"></a>\n",
    "\n",
    "Below we define a function that we will use to produce a color (RGB) image from a set of input images taken with three different filters. We'll use a tool from Astropy (`make_lupton_rgb`) 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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def create_rgb(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": [
    "#### 1.2.2. Parameters\n",
    "\n",
    "Let `afwDisplay` know that we want it to use `matplotlib` as our default display backend by setting the `setDefaultBackend()` function.\n",
    "\n",
    "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the parameters of `matplotlib.pyplot` to give us a large default size for an image, and set some other parameters to make the default plot style be color-blind friendly and otherwise make plots look nice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "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",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load the data to visualize\n",
    "\n",
    "To display an image, we must first load some data. The DP0.2 data set contains simulated images from the LSST DESC Data Challenge 2 (DC2) that have been reprocessed by the LSST Project using a more recent version of the LSST Science Pipelines. To access these data, we instantiate a Butler directing it to the `dp02` data repository configuration and the `2.2i/runs/DP0.2` collection.  For more information on the `Butler`, see [lsst.daf.butler](https://pipelines.lsst.io/modules/lsst.daf.butler/index.html) and subsequent tutorials in this series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "butler = Butler('dp02', collections='2.2i/runs/DP0.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a `Butler` instance generated, we can retrieve the desired calibrated exposure by telling the butler which visit and CCD (\"detector\") we wish to view. To do this, we define a dictionary with the required information. In this case, we access a single image from a specific visit (`192350`) and detector (`175`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dataId = {'visit': 192350, 'detector': 175, 'band': 'i'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "calexp = butler.get('calexp', **dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `calexp` that is returned by the Butler in the previous cell is an `ExposureF` Python object. Exposures are powerful representations of image data because they contain not only the image data, but also a variance image for uncertainty propagation, a bit mask image, and key-value metadata. In the next section, we will use AFWDisplay to visualize the image and mask associated with this Exposure. More documentation on accessing and visualizing an Exposure be found [here](https://pipelines.lsst.io/getting-started/display.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Basic image visualization\n",
    "\n",
    "### 3.1. Use AFWDisplay to visualize the image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now set to display the image. To do this, we:\n",
    "\n",
    "1. Create a `matplotlib.pyplot` figure using `plt.figure()` -- this will be familiar to anyone with experience using `matplotlib`.\n",
    "2. 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 `calexp.image` member of our calibrated image retrieved by the `butler`. We can then use `plt.show()` to display our figure.\n",
    "5. After we are done creating and displaying the image, we remove the underlying data from memory.\n",
    "\n",
    "All these tasks are best done within the same notebook cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(calexp.image)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 1: The `calexp` image, displayed in grayscale with a scale bar at right. Axes labels are in pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Use matplotlib imshow to visualize the image\n",
    "\n",
    "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` constructor 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`).\n",
    "\n",
    "To do this, we:\n",
    "1. Set the figure's projection to be the WCS of the `calexp`.\n",
    "2. Define the extent in pixel coordinates using the bounding box.\n",
    "3. Display the calexp image data array using the gray colormap (cmap) and use approximately the same min and max scale values as above.\n",
    "4. Add solid white grid lines.\n",
    "5. Label the axes, and show the plot.\n",
    "6. Remove the underlying data from memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.subplot(projection=WCS(calexp.getWcs().getFitsMetadata()))\n",
    "calexp_extent = (calexp.getBBox().beginX, calexp.getBBox().endX,\n",
    "                 calexp.getBBox().beginY, calexp.getBBox().endY)\n",
    "im = plt.imshow(calexp.image.array, cmap='gray', vmin=-200.0, vmax=400,\n",
    "                extent=calexp_extent, origin='lower')\n",
    "plt.grid(color='white', ls='solid')\n",
    "plt.xlabel('Right Ascension')\n",
    "plt.ylabel('Declination')\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 2: The same `calexp` image as in Figure 1, displayed with a different grayscale that increases the number of visible objects. Axes label and grid are in world coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.1. Use caution when interpreting imshow gridlines!\n",
    "\n",
    "In Figure 2, the more vertical white grid lines represent RA, and they are labled where they intersect the x-axis.\n",
    "In astronomy, the convention is to display images that are north-up.\n",
    "The rotation angle of the `calexp` displayed in Figure 2 happens to meet this convention, with east/west in RA along the x-axis, and north/south in Dec along the y-axis.\n",
    "\n",
    "However, the LSST Science Camera will obtain visits at a range of rotation angles.\n",
    "When visualizing a `calexp`, it is important to understand that it might be rotated such that east/west is along the y-axis instead.\n",
    "Here is a demonstraton of what that might look like.\n",
    "\n",
    "Define `dataId_2` and `calexp_2` for a visit that has a significant rotation angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId_2 = {'visit': 214463, 'detector': 171}\n",
    "calexp_2 = butler.get('calexp', **dataId_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the rotation angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(calexp_2.getInfo().getVisitInfo().boresightRotAngle.asDegrees())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the expected center of the image in RA, in hours-minutes-seconds format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbox = calexp_2.getBBox()\n",
    "wcs = calexp_2.getWcs()\n",
    "center = wcs.pixelToSky(bbox.centerX, bbox.centerY)\n",
    "c = SkyCoord(ra=center.getRa().asDegrees()*u.degree,\n",
    "             dec=center.getDec().asDegrees()*u.degree, frame='icrs')\n",
    "print(\"Rotated calexp's image center RA in h:m:s format is:\", c.ra.hms)\n",
    "del bbox, wcs, center, c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected image center in RA is thus 4h 8m 9s.\n",
    "\n",
    "Display the `calexp_2` with `matplotlib` `imshow` and overlay grid lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(projection=WCS(calexp_2.getWcs().getFitsMetadata()))\n",
    "calexp_extent = (calexp_2.getBBox().beginX, calexp_2.getBBox().endX,\n",
    "                 calexp_2.getBBox().beginY, calexp_2.getBBox().endY)\n",
    "ax.set_xlabel('RA (h:m:s)', fontsize=16)\n",
    "ax.set_ylabel('Dec (deg)', fontsize=16)\n",
    "ax.coords['ra'].set_ticks(number=12)\n",
    "ax.coords['dec'].set_ticks(number=12)\n",
    "im = plt.imshow(calexp_2.image.array, cmap='gray', vmin=-200.0, vmax=400,\n",
    "                extent=calexp_extent, origin='lower')\n",
    "plt.grid(color='white', ls='solid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 3: A visualization of rotated `calexp_2`, with grid lines. At first glance, the RA center of the image appears to be about 4h 8m 40s, which is far off the actual RA center of 4h 8m 9s obtained from the image WCS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**CAUTION:** Due to image rotation, the RA grid lines are actually the more horizontal lines in this display, but they are still labeled where they intersect the x-axis. Similarly, the Dec grid lines are the more vertical lines, yet still labeled where they intersect the y-axis. Contrast this with Figure 2, where the RA grid lines are more vertical (as expected for typicaly astronomy image display convention of north-up).\n",
    "\n",
    "Once the orientation of the image is understood, it becomes clearer that the RA center of the image is, indeed, about 4h 8m 9s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete `calexp_2` and the associated figure before moving on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "remove_figure(fig)\n",
    "del ax, im, calexp_extent\n",
    "del dataId_2, calexp_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. 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). \n",
    "One other component is the mask associated with the image. \n",
    "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": [
    "The `afwDisplay` package maps each bit in the mask plane to a specific display color. \n",
    "We can view this mapping using the code in the following cell. \n",
    "We can also use the `setMaskPlaneColor` method to change the colors that `afwDisplay` uses for each mask plane.\n",
    "\n",
    "Print the colors associated to each plane in the mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(\"Mask plane bit definitions:\\n\", display.getMaskPlaneColor())\n",
    "print(\"\\nMask plane methods:\\n\")\n",
    "help(display.setMaskPlaneColor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the image and mask plane side-by-side using matplotlib subplots.\n",
    "\n",
    "Use `plt.sca(ax[0])` to set the first axis as current, and then `plt.sca(ax[1])` to switch to the second axis.\n",
    "Using `plt.tight_layout()` with multi-axis figures helps to avoid axis overlap or excessive white spaces and results in a nicer-looking plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14, 7))\n",
    "plt.sca(ax[0])\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "display1.scale('linear', 'zscale')\n",
    "display1.mtv(calexp.image)\n",
    "plt.sca(ax[1])\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": [
    "> Figure 4: At left, the same `calexp` image as in Figures 1 and 2. At right, the mask plane for the `calexp` illustrates which pixels have a mask value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `afwDisplay` package also provides a nice interface for plotting the mask on top of the image using the `calexp.maskedImage`. \n",
    "The mask will also be plotted on top of the image if you pass the `calexp` itself to `mtv` (as is done later in this notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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": [
    "> Figure 5: The `calexp` image with the mask plane overlaid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To investigate the mask in a bit more detail, follow the same steps as above to display the image, but add a few modifications:\n",
    "\n",
    "1. Explicitly set the transparency of the overplotted mask to 10% (0 = opaque, 100 = transparent).\n",
    "2. Change the color of the 'DETECTED' mask plane from 'blue' (as above) to 'green' (i.e., all pixels associated with detected objects).\n",
    "3. Pass the full `calexp` object to `mtv` instead of just the image plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.setMaskTransparency(10)\n",
    "display.setMaskPlaneColor('DETECTED', 'green')\n",
    "display.mtv(calexp)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 6: Similar to Figure 4, but with the mask transparency reduced to 10% (mostly opaque) and the color representing 'DETECTED' pixels changed from blue to green."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4. More information about `lsst.afw.display`\n",
    "\n",
    "To find more information about `lsst.afw.display`, print the method list to see what's available. \n",
    "The next cell will print `lsst.afw.display` methods to the screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# help(display.scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# help(display.mtv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Congratulations!** We've plotted an image in various ways using `lsst.afw.display`!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Extract cutout images\n",
    "\n",
    "When a closer look at images in the vicinity of a star or galaxy is needed, visualizing a small cutout instead of the full image is more effective.\n",
    "\n",
    "### 4.1. A cutout from a `calexp`\n",
    "\n",
    "In a science use case, the coordinates of the object of interest may be known already, but for this example choose pixels near an extended source in the image above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wcs = calexp.getWcs()\n",
    "radec = wcs.pixelToSky(2250, 700)\n",
    "ra, dec = radec.getRa().asDegrees(), radec.getDec().asDegrees()\n",
    "print(ra, dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the same visit and detector as were used for the `dataId` in Section 3.\n",
    "\n",
    "> **Notice:** In order to make a cutout of a `calexp`, the RA, Dec, visit, *and* detector must all be specified (recall that visit uniquely identifies band, so band does not need to be explicitly identified).\n",
    "This is different from making a cutout with a `deepCoadd` (Section 4.2 below), which only requires RA, Dec, and band.\n",
    "This is because there are multiple `calexp` images for a given coordinate and band, but only one `deepCoadd` per band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "visit = 192350\n",
    "detector = 175"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `cutout_calexp` function defined in Section 1 to extract a 301x301 pixel cutout.\n",
    "\n",
    "> **Optional:** print the help documentation for the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# help(cutout_calexp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_cutout_calexp = cutout_calexp(butler, ra, dec, visit, detector, cutoutSideLength=301)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the cutout from the `calexp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(my_cutout_calexp.image)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 7: A cutout from a `calexp` image, displayed in grayscale, with scale bar at right and axes in pixels. The three brightest objects are galaxies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del wcs, radec, ra, dec, visit, detector, my_cutout_calexp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. A cutout from a `deepCoadd` image\n",
    "\n",
    "Up until now, we have been examining an image from a single CCD detector taken during a single visit. For the rest of this notebook, we will switch to examining coadded images made up of multiple exposures that have been combined. \n",
    "Let's start by taking a look at what a full 4k x 4k pixel coadd \"patch\" image looks like. \n",
    "\n",
    "We start by grabbing a `deepCoadd` image and displaying it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dataId = {'tract': 4431, 'patch': 17, 'band': 'i'}\n",
    "datasetType = 'deepCoadd'\n",
    "coadd = butler.get(datasetType, **dataId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(coadd.image)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 8: A `deepCoadd` patch image, displayed in grayscale with scale bar at right, showing a rich galaxy cluster in the lower-left corner."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image above is displaying pixel coordinates (note that the coadd patch is part of a larger coadd image called a \"tract\", so the pixel values do not start at 0,0), but in general it is more useful to be able to select a region based on RA, Dec coordinates. To do this, we'll use the world coordinate system (WCS) object associated with the image.\n",
    "\n",
    "Extract the WCS solution, which provides the mapping between XY pixel values and sky coordinates, and print it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wcs = coadd.getWcs()\n",
    "print(wcs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cluster seems to be centered at about (X, Y) ~ (12500, 8500).\n",
    "Use the \"pixelToSky\" method of the WCS to get the sky coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "radec = wcs.pixelToSky(12500, 8500)\n",
    "ra, dec = radec.getRa().asDegrees(), radec.getDec().asDegrees()\n",
    "print(ra, dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the RA, Dec coordinates of the cluster, we would like to grab a small cutout of the coadded image at this location. To do this, we've defined a user function [`cutout_coadd` at the beginning of the notebook](#cutout_coadd). This function extracts a cutout from a deep coadd image at a given RA, Dec position and desired image size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "help(cutout_coadd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call this function to extract a cutout image that is centered on the galaxy cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cutout_image = cutout_coadd(butler, ra, dec, band='i', 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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout_image.image)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 9: A cutout showing a \"zoom in\" on Figure 7."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Overplotting catalog sources\n",
    "\n",
    "Now that we've been able to examine the coadd image and make a cutout, we may be interested in the corresponding catalog of astronomical sources that was extracted from this image by the LSST Science Pipelines. \n",
    "The TAP service is the recommended way to retrieve catalog data for a notebook, and there are examples on how to do this in previous [tutorials](https://github.com/rubin-dp0/tutorial-notebooks).\n",
    "\n",
    "This section demonstrates how to use TAP to retrieve the catalog data within a polygon defined by the corners of the coadd image.\n",
    "\n",
    "The following cell extracts the XY values of the corners of the cutout image and converts them to RA, Dec, which are then used as spatial constraints in the query below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wcs = cutout_image.getWcs()\n",
    "\n",
    "x0 = float(cutout_image.getX0())\n",
    "y0 = float(cutout_image.getY0())\n",
    "width = cutout_image.getWidth()\n",
    "height = cutout_image.getHeight()\n",
    "\n",
    "xcorners = [x0, x0+width, x0+width, x0]\n",
    "ycorners = [y0, y0, y0+width, y0+width]\n",
    "\n",
    "ra_corners = []\n",
    "dec_corners = []\n",
    "\n",
    "for i in range(len(xcorners)):\n",
    "    radec = wcs.pixelToSky(xcorners[i], ycorners[i])\n",
    "    ra_corners.append(radec.getRa().asDegrees())\n",
    "    dec_corners.append(radec.getDec().asDegrees())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instantiate the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a query statement for `Objects` within the bounding box corners.\n",
    "\n",
    "> **Notice:** You may be familiar with the `BETWEEN ... AND ...` format for ADQL spatial queries, and tempted to use it here to select on RA and Dec ranges. However, the Qserv database does not know how to deal with this construction efficiently, so you are _much_ better off constructing a polygon from the corners of the image as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, coord_ra, coord_dec, x, y, tract, patch \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \" +\\\n",
    "        \"POLYGON('ICRS', \" + str(ra_corners[0]) + \", \" + str(dec_corners[0]) + \", \" +\\\n",
    "        str(ra_corners[1]) + \", \" + str(dec_corners[1]) + \", \" +\\\n",
    "        str(ra_corners[2]) + \", \" + str(dec_corners[2]) + \", \" +\\\n",
    "        str(ra_corners[3]) + \", \" + str(dec_corners[3]) + \")) = 1 AND \" +\\\n",
    "        \"tract = \" + str(dataId['tract']) + \" AND patch = \" + str(dataId['patch'])\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the query as an asynchronous job, then fetch the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the results to an Astropy table, then take a look at the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tab = results.to_table()\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the image cutout, and overplot the sources on the image. Use display buffering to avoid re-drawing the image after each source is plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout_image.image)\n",
    "with display.Buffering():\n",
    "    for i in range(len(tab)):\n",
    "        display.dot('+', tab[i]['x'], tab[i]['y'], ctype=afwDisplay.RED)\n",
    "        display.dot('o', tab[i]['x'], tab[i]['y'], size=20, ctype='orange')\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 10: The same `deepCoadd` cutout as in Figure 8, but with red crosses and orange circles marking the location of `Objects`.\n",
    "\n",
    "Note that since the object catalog is assembled from the multi-band coadd images, not every source will necessarily be detected in our i-band image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Plot an RGB composite image\n",
    "\n",
    "We've found a pretty galaxy cluster, but what if we 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 and then assign those images to the RGB channels of a color image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the `MultiBandExposure` object type from the LSST Science Pipelines `afw.image` package. \n",
    "However, one could also create a similar function that simply accepts three separate images without combining them into a MultiBand object.\n",
    "Notice that the `MultibandExposure.fromExposures` function needs a list of images and filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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",
    "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). We will use the function [`create_rgb`](#create_rgb) that we defined at the [top of this notebook](#create_rgb).\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 that better highlights certain features. You can experiment with the \"scale = []\" keyword, or including different bands besides _gri_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(20, 20), nrows=1, ncols=2)\n",
    "\n",
    "rgb_original = create_rgb(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",
    "ax[1].set_title('re-scaled', fontsize=30)\n",
    "rgb_scaled = create_rgb(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",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Figure 11: Color images built from slightly larger cutouts that are centered on the same coordinates as in Figures 8 and 9. At left, the default RGB scaling is used; at right, a scaling that gives more weight to blue is used, resulting in a cooler palette. "
   ]
  },
  {
   "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": [
    "## 7. Exercises for the learner\n",
    "\n",
    "1. Create `coadd` and `calexp` cutouts for the same RA, Dec, and band, and compare them side-by-side.\n",
    "2. Write your own function to create a cutout of a color image.\n",
    "\n",
    "## 8. 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.9"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}