IAUGA_2024_images.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br><b>Image display demo for the IAU General Assembly LSST RSP workshop, August 2024</b> <br>\n",
    "Rubin Science Platform deployment: data.lsst.cloud <br>\n",
    "LSST Science Pipelines version: Weekly 2024_16 <br>\n",
    "Last verified to run: Tue Aug 13 2024 <br>\n",
    "Contact author: Christina Williams <br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "Learn more about the Data Preview 0 (DP0) data set (e.g., image types, catalog tables)\n",
    "and the LSST Science Pipelines and Rubin Science Platform functionality using the\n",
    "documentation and resources available at <a href=\"https://dp0.lsst.io/\">dp0.lsst.io</a>.\n",
    "\n",
    "<br>\n"
   ]
  },
  {
   "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 extensions and metadata, and create cutouts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**LSST Data Products:** Images (calexp and deepCoadd image types), table that records the location and the type of image (ivoa.ObsCore table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**LSST Python Packages:** lsst.afw.display, lsst.daf.butler, lsst.geom, lsst.afw.image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial is based on multiple tutorials developed by many people, including tutorial notebooks DP02_01, DP02_02a,b,c, and DP02_03a,b that are published by the Rubin Community Science Team [at the DP0.2 Tutorial-Notebooks github repository](https://pipelines.lsst.io/getting-started/display.html). Rendered versions (including output) are available [here](https://dp0-2.lsst.io/tutorials-examples/rendered-tutorial-notebooks.html).\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 introduce users to accessing imaging data products produced by the LSST pipelines, their format, and basic image display and manipulation. \n",
    "\n",
    "It also introduces LSST pipeline packages that are useful for basic access:  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 [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 display the image data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### 1.1. Package imports\n",
    "\n",
    "Below, the [`matplotlib.pyplot`](https://matplotlib.org/) library 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 `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": [
    "# general-use python packages\n",
    "import matplotlib.pyplot as plt\n",
    "import gc\n",
    "\n",
    "# Astropy packages\n",
    "from astropy import units as u\n",
    "from astropy.wcs import WCS\n",
    "\n",
    "# LSST Pipeline packages for data access:\n",
    "from lsst.daf.butler import Butler\n",
    "from lsst.rsp import get_tap_service\n",
    "\n",
    "# LSST Pipeline package for dealing with sky coordinates\n",
    "import lsst.geom as geom\n",
    "\n",
    "# LSST Pipieline packages for image manipulation and display\n",
    "import lsst.afw.display as afwDisplay\n",
    "from lsst.afw.image import ExposureF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2.1. Functions\n",
    "\n",
    "After producing a plot, the matplotlib package will store the data array associated with an image in memory. 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": [
    "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. Intro to Image Data\n",
    "\n",
    "This tutorial focuses on two common types of images people will interact with: `calexp` and `deepCoadd`.\n",
    "\n",
    "<b>calexp: A calibrated exposure (single image in a single filter).</b>\n",
    "\n",
    "Calexp files are uniquely identified by `visit` (when observation was taken) and `detector` in the camera.  \n",
    "\n",
    "<img src=\"figures/img_demonstration.png\" width=\"400\">\n",
    "Figure 1: Example of many calexp footprints (magenta) that are combined into deepCoadds (gray)\n",
    "<br><br>\n",
    "<b>deepCoadd: A combination of single images into a deep stack or coadd.</b>\n",
    "\n",
    "The LSST Science Pipelines processes and stores `deepCoadd` images in tracts and patches.\n",
    "\n",
    "tract: A portion of sky within the LSST all-sky tessellation (sky map); divided into patches.\n",
    "\n",
    "patch: A quadrilateral sub-region of a tract, of a size that fits easily into memory on desktop computers.\n",
    "\n",
    "<img src=\"figures/dp02_patch_tract.png\" width=\"700\">\n",
    "Figure 2: Example tract divided into patches in DP0.2 (Abolfathi et al. 2021)\n",
    "<br><br>\n",
    "\n",
    "To retrieve and display an image at a desired coordinate, users have to specify their image type, and additional parameters, which make up whats called a `dataId` : tract, patch and filter are needed in the case of `deepCoadd` ; only the visit and detector are needed in the case of `calexp`.\n",
    "\n",
    "## 2.1. Retrieve an image using the LSST Butler\n",
    "\n",
    "The Butler (<a href=\"https://pipelines.lsst.io/modules/lsst.daf.butler/index.html\">documentation</a>) is an LSST Science Pipelines software package to fetch LSST data without having to know its location or format. This is the preferred method for image retrieval using the RSP Notebook aspect.\n",
    "\n",
    "The DP0.2 data set contains simulated images from the LSST Dark Energy Science Collaboration (DESC) Data Challenge 2 (DC2) that have been reprocessed by the LSST Project using a more recent version of the LSST Science Pipelines.  \n",
    "\n",
    "### 2.1.1 Create an instance of the butler \n",
    "Specify the set of imaging to use in the RSP. Request simulated DP0.2 imaging by directing the Butler to the `dp02` data repository configuration, and the collection `2.2i/runs/DP0.2` specifies which processing run using the LSST pipelines to use. \n",
    "\n",
    "It will return an informative statement about credentials being found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butler = Butler('dp02', collections='2.2i/runs/DP0.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.2 Identify and retrieve a deepCoadd\n",
    "There is a nice galaxy cluster in the DP0.2 simulated data at RA = 03h42m59.0s, Dec = -32d16m09s (in degrees, 55.745834, -32.269167).\n",
    "\n",
    "Use the LSST `geom` package to define a `SpherePoint` (an LSST pipeline object to hold coordinates) for the cluster's coordinates (<a href=\"https://pipelines.lsst.io/modules/lsst.geom/index.html\">lsst.geom documentation</a>).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_ra = 55.745834\n",
    "my_dec = -32.269167\n",
    "\n",
    "my_spherePoint = geom.SpherePoint(my_ra*geom.degrees,\n",
    "                                  my_dec*geom.degrees)\n",
    "print(my_spherePoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the Butler to retrieve the DP0.2 sky map. The sky map is simply a reference of how the all sky area is divided up into patch and tracts (see Figure 2). Use it to identify the tract and patch for the cluster's coordinates (<a href=\"http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/skymap.html\">skymap documentation</a>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skymap = butler.get('skyMap')\n",
    "\n",
    "tract = skymap.findTract(my_spherePoint)\n",
    "patch = tract.findPatch(my_spherePoint)\n",
    "\n",
    "my_tract = tract.tract_id\n",
    "my_patch = patch.getSequentialIndex()\n",
    "\n",
    "print('my_tract: ', my_tract)\n",
    "print('my_patch: ', my_patch)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the Butler to retrieve the i-band `deepCoadd` for the tract and patch. To do this, define a dictionary called `dataId` with the required information (for deepCoadds, to uniquely identify the image, the user must provide 3 things: the filter (band), patch, and tract). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = {'band': 'i', 'tract': my_tract, 'patch': my_patch}\n",
    "my_deepCoadd = butler.get('deepCoadd', dataId=dataId)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.3 LSST image format\n",
    "\n",
    "The `my_deepCoadd` returned by the Butler contains the science image, and also has a number of extensions containing other data or information about the image. This includes the science image, a mask plane containing the masked pixels and why they were masked, the variance image, and spatially varying PSF information.\n",
    "\n",
    "These examples will focus on just the science and mask extensions (see the [DP0.2 Tutorial-Notebooks github repository](https://pipelines.lsst.io/getting-started/display.html) for additional demonstrations). \n",
    "\n",
    "The extension containing the science image and masks can be returned using the `image` and `mask` attributes. One can retrieve just an array containing the pixels values using `.array` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deepCoadd_image = my_deepCoadd.image\n",
    "deepCoadd_pixels = my_deepCoadd.image.array\n",
    "deepCoadd_mask = my_deepCoadd.mask\n",
    "deepCoadd_var = my_deepCoadd.variance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 3. Basic image visualization\n",
    "\n",
    "### 3.1 Use AFWDisplay to visualize the image\n",
    "Image data retrieved with the Butler is an `ExposureF` (F=float) Python object. `ExposureF` is the native format for imaging data taken with Rubin and processed by LSST Pipelines. They 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 object be found [here](https://pipelines.lsst.io/getting-started/display.html).\n",
    "\n",
    "The image contained in the `ExposureF` object can be displayed several different ways. A simple option is to use the LSST Science Pipelines package `afwDisplay`. There is some <a href=\"https://pipelines.lsst.io/modules/lsst.afw.display/index.html\">documentation for afwDisplay</a> available, and [other DP0.2 Tutorial-Notebooks](https://pipelines.lsst.io/getting-started/display.html) (DP02_03a,b)  go into more detail about all the display options.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.1 Interactive Display using afwDisplay\n",
    "\n",
    "First, explore the image using an interactive display called `firefly`, which is the RSP's interface for interactive image display and manipulation. \n",
    "\n",
    "1. 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 `firefly_display`.\n",
    "2. Display the image by providing the `mtv()` method the `deepCoadd.image` member of our calibrated image retrieved by the Butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('firefly')\n",
    "firefly_display = afwDisplay.Display(frame=1)\n",
    "firefly_display.mtv(deepCoadd_image)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.2 Static Image Display using afwDisplay\n",
    "\n",
    "It is also possible to plot statically in the notebook. To do this, set the backend of `afwDisplay` to `matplotlib`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The following code cell creates a `matplotlib.pyplot` figure; aliases `lsst.afw.display.Display` as `display`; \n",
    "sets the image stretch for the pixel shading (set by the algorithm asinh -- familiar from SDSS images -- with a range of values set by zscale); displays the image data using the `mtv` method; and turns on the x and y axes labels (pixel coordinates).\n",
    "\n",
    "At the end, remove the figure data to free some memory.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "display = afwDisplay.Display(1)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(my_deepCoadd.image)\n",
    "plt.gca().axis('on')\n",
    "plt.show()\n",
    "\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.3 More information and help\n",
    "\n",
    "To find more information, the next cell will print the `lsst.afw.display` methods that are available for this display type 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": [
    "To learn more about the afwDisplay package and its tasks (e.g. `mtv`), use the help function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(display.mtv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 3.2 Working with World Coordinate Systems (WCS)\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 `deepCoadd`.\n",
    "2. Define the extent in pixel coordinates using the bounding box.\n",
    "3. Display the deepCoadd image data array using the gray colormap (cmap)\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",
    "\n",
    "my_deepCoadd_WCS = WCS(my_deepCoadd.getWcs().getFitsMetadata())\n",
    "\n",
    "plt.subplot(projection = my_deepCoadd_WCS)\n",
    "\n",
    "my_deepCoadd_extent = (my_deepCoadd.getBBox().beginX, my_deepCoadd.getBBox().endX,\n",
    "                       my_deepCoadd.getBBox().beginY, my_deepCoadd.getBBox().endY)\n",
    "\n",
    "im = plt.imshow(my_deepCoadd.image.array, cmap='gray', vmin=-2.0, vmax=4,\n",
    "                extent=my_deepCoadd_extent, origin='lower')\n",
    "\n",
    "plt.grid(color='white', ls='solid')\n",
    "plt.xlabel('Right Ascension')\n",
    "plt.ylabel('Declination')\n",
    "plt.show()\n",
    "\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Use afwDisplay to visualize the image and mask plane\n",
    "\n",
    "The `deepCoadd` returned by the Butler contains more than just the image pixel values (see this [image 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": [
    "Now plot the image and mask side-by-side using matplotlib.\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",
    "\n",
    "# set active plot to left panel\n",
    "plt.sca(ax[0])\n",
    "\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "display1.scale('linear', 'zscale')\n",
    "display1.mtv(my_deepCoadd.image)\n",
    "\n",
    "# right panel\n",
    "plt.sca(ax[1])\n",
    "\n",
    "display2 = afwDisplay.Display(frame=fig)\n",
    "display2.mtv(my_deepCoadd.mask)\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `afwDisplay` package also provides a nice interface for plotting the mask on top of the image using the `my_deepCoadd.maskedImage`. \n",
    "The mask will also be plotted on top of the image if you pass the `my_deepCoadd` itself to `mtv` (as is done later in this notebook). The below cell will print the types of pixel masks that are available, and plot one mask rather than all masked pixels (as done above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "\n",
    "for maskName, maskBit in my_deepCoadd.mask.getMaskPlaneDict().items(): \n",
    "    print(f'{maskName:18s} (bit=2**{maskBit:02d})')\n",
    "    display.setMaskTransparency(100,maskName)\n",
    "\n",
    "display.setMaskTransparency(10,'DETECTED') \n",
    "display.setMaskPlaneColor('DETECTED', 'blue')\n",
    "\n",
    "# optionally plot any pixel that was at the image edge of a calexp\n",
    "#display.setMaskTransparency(10,'SENSOR_EDGE') \n",
    "#display.setMaskPlaneColor('SENSOR_EDGE', 'yellow')\n",
    "\n",
    "display.scale('linear', 'zscale')\n",
    "display.mtv(my_deepCoadd.maskedImage)\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To investigate the mask in a bit more detail, 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 `my_deepCoadd` object to `mtv` instead of just the image plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.setMaskTransparency(80)\n",
    "display.setMaskPlaneColor('DETECTED', 'green')\n",
    "display.mtv(my_deepCoadd)\n",
    "plt.show()\n",
    "\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "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 = my_deepCoadd.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",
    "cluster_ra, cluster_dec = radec.getRa().asDegrees(), radec.getDec().asDegrees()\n",
    "print(cluster_ra, cluster_dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 4. Accessing Images via TAP query: The ObsCore table\n",
    "\n",
    "It is also possible to use TAP queries to access imaging information from an LSST catalog, and does not require the use of LSST Pipelines like the Butler does. A use case for this would be if one wants to access images from outside the RSP.\n",
    "\n",
    "The LSST `ObsCore` table contains information about the images\n",
    "stored in the LSST's [data butler](https://pipelines.lsst.io/modules/lsst.daf.butler/index.html) registry (for more information, see data management document, [DMTN-236](https://dmtn-236.lsst.io/)). The IVOA prefix in the catalog's title `ivoa.ObsCore` stands for International Virtual Observatory Alliance, which sets standards for data storage (the Virtual Observatory; VO) and Rubin's data products are consistent with these VO standards.\n",
    "\n",
    "For this part of the tutorial, import some standard VO packages that will enable data access via VO standards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PyVO packages\n",
    "import pyvo\n",
    "from pyvo.dal.adhoc import DatalinkResults, SodaQuery\n",
    "\n",
    "# more LSST Pipeline packages\n",
    "from lsst.resources import ResourcePath\n",
    "from lsst.afw.fits import MemFileManager\n",
    "from lsst.afw.image import ExposureF\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like Tutorial 2 on catalogs, we initialize the TAP service to allow queries to the `ivoa.ObsCore` table. The second line will store an authentication that is returned by TAP that allows data access. Authentication is needed to access the imaging (if proprietary) from outside the RSP environment (for more information see https://rsp.lsst.io/guides/auth/creating-user-tokens.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")\n",
    "auth_session = service._session"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-08-03T23:59:19.228126Z",
     "iopub.status.busy": "2024-08-03T23:59:19.227688Z",
     "iopub.status.idle": "2024-08-03T23:59:19.245715Z",
     "shell.execute_reply": "2024-08-03T23:59:19.244228Z",
     "shell.execute_reply.started": "2024-08-03T23:59:19.228074Z"
    }
   },
   "source": [
    "Retrieve the names of the columns and their data types, descriptions and units,\n",
    "and display these properties as a pandas table (allows display of all rows/columns). The output is a list of the image metadata that is available in the `ivoa.ObsCore` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT column_name, datatype, description, unit \" \\\n",
    "        \"FROM tap_schema.columns \" \\\n",
    "        \"WHERE table_name = 'ivoa.ObsCore'\"\n",
    "\n",
    "results = service.search(query)\n",
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Useful metadata includes `s_region`, which is the bounding box (spatial extent) of the image. This can be used to identify which images cover the ra/dec of interest.\n",
    "\n",
    "The following query selects all columns from the `ivoa.ObsCore` table\n",
    "for all images that contain the coordinate for the galaxy cluster within the `s_region`. (In other words, where the statment that the `s_region` contains \n",
    "the ra & dec of the point is `True` or `= 1`.)\n",
    "\n",
    "For large queries like this one, we submit \"asynchronously\" and wait for the job to complete before continuing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT * FROM ivoa.ObsCore \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', \" + str(my_ra) + \\\n",
    "        \", \" + str(my_dec) + \"), s_region) = 1\"\n",
    "\n",
    "# Define the service job using the query, then run it.\n",
    "job = service.submit_job(query)\n",
    "job.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait until the job is completed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the job phase is \"COMPLETED\", fetch (retrieve) the results,\n",
    "and print the number of rows returned.\n",
    "This is the number of images found that overlap the specified point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results), \"images covering this ra/dec are found\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Retrive image data (headers and pixel data), for imaging. This time, we will retrieve a `calexp` image (which are single visit exposures on sky; see Section 2). \n",
    "\n",
    "`calexp` (calibrated exposures) are classified as calibration level (`calib_level`) of 2. `calib_level` is an IVOA convention to classify the different level of processing that has been applied to the observed dataset. \n",
    "\n",
    "\n",
    "In the RSP, the dataset calib_level numbers correspond to:\n",
    "\n",
    "    1. Raw images (not science-ready)\n",
    "    2. Calibrated single-visit exposures (calexp)\n",
    "    3. Coadded images and difference images (e.g. deepCoadd)\n",
    "\n",
    "\n",
    "For Rubin/LSST, we will down-sample the query to only return `calexp` images (i.e. `calib_level`=2) taken in the i-band. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT * FROM ivoa.ObsCore \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', \" + \\\n",
    "        str(my_ra) + \", \" + str(my_dec) + \"), s_region) = 1 \"\\\n",
    "        \"AND lsst_band = 'i' \"\\\n",
    "        \"AND calib_level = 2\"\n",
    "\n",
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "results = job.fetch_result()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Option to view results as an `astropy` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ra/dec is covered by 151 `calexp` images.  For this demonstration only one of the 151 images is needed. So, extract any row into `results_0` (in this case we will use row 100)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_0 = results[100]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print selected metadata:\n",
    " * dataproduct subtype\n",
    " * CCD visit identifier (`lsst_ccdvisitid`)\n",
    " * visit identifier and detector number\n",
    " * central RA and Dec\n",
    " * band (filter)\n",
    " * exposure start time `t_min` (a modified julian date, MJD)\n",
    " * `s_region` (defines the corner coordinates of the image; i.e. its extent)\n",
    " * a unique `obs_id` in the TAP service"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_0['dataproduct_subtype'])\n",
    "print(results_0['lsst_ccdvisitid'])\n",
    "print(\"visit number = \", results_0['lsst_visit'], \n",
    "      \"detector number = \", results_0['lsst_detector'])\n",
    "print('%7.4f %8.4f' % (results_0['s_ra'], results_0['s_dec']))\n",
    "print(results_0['lsst_band'])\n",
    "print(results_0['t_min'])\n",
    "print(results_0['s_region'])\n",
    "print(results_0['obs_id'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1 Access the image via a Datalink URL\n",
    "\n",
    "Data is accessible via an `access_url` and the format of the data is stored in `access_format`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_0['access_url'])\n",
    "print(results_0['access_format'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, the format indicates that the image information is stored as a `votable` which is accessible via Datalink.\n",
    "\n",
    "Obtain the table `dl_results` using Datalink and the `auth_session` that was defined earlier to authenticate data access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dl_results = DatalinkResults.from_result_url(results_0['access_url'],\n",
    "                                             session=auth_session)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Option to view `dl_results` as an astropy table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dl_results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For security reasons the `dl_results` table holds a *different* `access_url` that\n",
    "provides direct access to the image - define it as `image_url` and print it to see how it is different from the `access_url` above. We need the datalink URL, which was provided once authenticated, to access the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_url = dl_results['access_url'][0]\n",
    "print(image_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the LSST pipelines are still in development, there is not a way to read the `ExposureF` directly from a url currently. To access the image: first do the extra step of allocating memory for the image data (stored in `mem`) before displaying the image.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fits_data = ResourcePath(image_url).read()\n",
    "\n",
    "mem = MemFileManager(len(fits_data))\n",
    "mem.setData(fits_data, len(fits_data))\n",
    "exposure = ExposureF(mem)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(exposure.image)\n",
    "plt.show()\n",
    "\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Optional: Creating cutouts of LSST images using the cutout tool\n",
    "\n",
    "Image cutouts are an important tool that enables a faster way to visually inspect science targets. Below we demonstrate the current way to perform image cutouts, using the IVOA-compatible cutout tool.\n",
    "\n",
    "These pipelines are still under development, in the future its likely these will be packaged into more user-fiendly functions. We may develop new ways that would be recommended ways for users to do these operations.\n",
    "\n",
    "(Clarification note: the `dataproduct_subtype` called `deepCoadd_calexp` is similar to `deepCoadd` but with a small background removed to distinguish it from `deepCoadd`. During operations with real data, deepCoadds will be called `deepCoadd_calexp`. This is confusing nomenclature for DP0.2, but will be clarified in the future)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT access_format, access_url, dataproduct_subtype, \" + \\\n",
    "    \"lsst_patch, lsst_tract, lsst_band, s_ra, s_dec  \" + \\\n",
    "    \"FROM ivoa.ObsCore WHERE dataproduct_type = 'image' \" + \\\n",
    "    \"AND obs_collection = 'LSST.DP02' \" + \\\n",
    "    \"AND dataproduct_subtype = 'lsst.deepCoadd_calexp' \" + \\\n",
    "    \"AND lsst_tract = \" + str(my_tract) + \" \" + \\\n",
    "    \"AND lsst_patch = \" + str(my_patch) + \" \" + \\\n",
    "    \"AND lsst_band = 'i' \"\n",
    "\n",
    "results = service.search(query)\n",
    "results.to_table().show_in_notebook()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# obtain the access_url for the data in the LSST database, \n",
    "# as we did for the calexp \n",
    "access_url = results[0].getdataurl()\n",
    "\n",
    "# use the access_url to get the datalink url pointing directly to data\n",
    "dl_results = DatalinkResults.from_result_url(access_url,\n",
    "                                             session=auth_session)\n",
    "\n",
    "f\"Datalink status: {dl_results.status}. Datalink service url: {access_url}\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below calls a cutout service (it is a function enabled by the `pyvo` package as part of a virtual observatory). The cutout shape is defined by the `circle` function, which will create a small square image with edge size = diameter of the defined circle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sq = SodaQuery.from_resource(dl_results,\n",
    "                             dl_results.get_adhocservice_by_id(\"cutout-sync\"),\n",
    "                             session=auth_session)\n",
    "\n",
    "# build the shape of the cutout. Circle returns a square cutout with \n",
    "# edge size = diameter of the subtended circle\n",
    "sphereRadius = 0.05 * u.deg\n",
    "\n",
    "sq.circle = (my_spherePoint.getRa().asDegrees() * u.deg,\n",
    "             my_spherePoint.getDec().asDegrees() * u.deg,\n",
    "             sphereRadius)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like earlier, store the resulting cutout in memory (`mem`) and then plot the returned cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# store the cutout data in memory:\n",
    "cutout_bytes = sq.execute_stream().read()\n",
    "mem = MemFileManager(len(cutout_bytes))\n",
    "mem.setData(cutout_bytes, len(cutout_bytes))\n",
    "exposure = ExposureF(mem)\n",
    "\n",
    "# plot the image:\n",
    "fig, ax = plt.subplots()\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(exposure.image)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## 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.11.7"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}