DP02_13a_Image_Cutout_SciDemo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# Using the image cutout tool with DP0.2\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\" alt=\"Rubin Observatory logo, a graphical representation of turning stars into data.\">\n",
    "<br>\n",
    "Contact author: Christina Williams <br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Piplines version: Weekly 2024_50 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** This notebook demonstrates how to use the Rubin Image Cutout Service."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Run the Rubin Image Cutout Service for visual inspection of small cutouts of LSST images."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** Images (deepCoadd, calexp), catalogs (objectTable, diaObject, truthTables, ivoa.ObsCore)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** PyVO, lsst.rsp.get_tap_service, lsst.pipe.tasks.registerImage, lsst.afw.display"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:**\n",
    "This notebook builds on technical exploration of the cutout service developed by Leanne Guy (in preparation), and builds on an earlier notebook written by Alex Drlica-Wagner <a href=\"https://github.com/rubin-dp0/cst-dev/blob/main/ADW_gravelpit/TAP_Image_Access.ipynb\">linked here</a>. This notebook additionally includes contributed functions by Jeff Carlin and Ryan Lau for time-variable sources. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28e91cbf-ab7f-4e26-9276-b00299d6065e",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. 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",
   "id": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook will teach how to use the cutout service (which enables users to retrieve small image cutouts from LSST data) by demonstration of a few basic science applications. Since LSST images are quite large, in some cases it is desirable to be able to perform image cutouts on the server side to avoid transferring large amounts of data. This can speed up visual inspection and other analyses requiring images of individual objects or fields that may be much smaller than the patch and tract sizes that are created by the LSST camera and pipelines.\n",
    "\n",
    "The <a href=\"https://www.ivoa.net/\">International Virtual Observatory Alliance (IVOA)</a> co-ordinates the community efforts of astronomical missions and archives to develop and maintain the Virtual Observatory (VO) standards. The VO standards enable interoperability between astronomical archives. IVOA also provides a <a href=\"https://www.ivoa.net/documents/SODA/20170517/REC-SODA-1.0.html\">Server-side Operations for Data Access (SODA)</a>, which is a protocol for remote data processing operations. This protocol allows users to perform computations (pixel operations, image transformations, etc) on the remote server, which avoids unnecssary data movement. The LSST architecture has a \"VO-first\" approach, meaning that VO standards are implemented in all applicable services, enabling the use of VO tools such as the image cutout service to access LSST data. The procedure is to identify the remote web location of the image of interest (called a datalink), and use a web service that creates a cutout from that linked data remotely, to transfer and save locally on the Rubin Science Platform.\n",
    "\n",
    "Further details and information can be found at <a href=\"https://www.ivoa.net/documents/DataLink/20211115/WD-DataLink-1.1-20211115.html\"> the IVOA data link documentation</a>, where it says Access Data Services. Rubin-specific documentation for these can also be found in <a href=\"https://dmtn-238.lsst.io/\">this document describing the RSP DataLink service implementation strategy</a>. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1 Package Imports\n",
    "\n",
    "Import general python packages, LSST packages, PyVO packages, and Astropy packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import uuid\n",
    "import os\n",
    "import glob\n",
    "import math\n",
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "from IPython.display import Image as dimg\n",
    "\n",
    "import lsst.geom as geom\n",
    "import lsst.resources\n",
    "import lsst.afw.display as afwDisplay\n",
    "from lsst.afw.image import Exposure, ExposureF\n",
    "from lsst.pipe.tasks.registerImage import RegisterConfig, RegisterTask\n",
    "from lsst.rsp import get_tap_service\n",
    "from lsst.rsp.utils import get_access_token\n",
    "from lsst.afw.fits import MemFileManager\n",
    "\n",
    "import pyvo\n",
    "from pyvo.dal.adhoc import DatalinkResults, SodaQuery\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters\n",
    "\n",
    "Set the backend for `afwDisplay` to `matplotlib`.\n",
    "Set parameters for plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "552c4cac-feca-4d51-a92f-87d5779a3300",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "\n",
    "params = {'axes.labelsize': 18,\n",
    "          'font.size': 18,\n",
    "          'legend.fontsize': 12,\n",
    "          'xtick.major.width': 2,\n",
    "          'xtick.minor.width': 1,\n",
    "          'xtick.major.size': 10,\n",
    "          'xtick.minor.size': 4,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 2,\n",
    "          'axes.linewidth': 2,\n",
    "          'axes.labelweight': 2,\n",
    "          'axes.titleweight': 2,\n",
    "          'ytick.major.width': 2,\n",
    "          'ytick.minor.width': 1,\n",
    "          'ytick.major.size': 10,\n",
    "          'ytick.minor.size': 4,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [6, 6],\n",
    "          'figure.facecolor': 'White'\n",
    "          }\n",
    "plt.style.use('tableau-colorblind10')\n",
    "\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "522f9661-2534-418a-bea4-5e6bbf5eb615",
   "metadata": {},
   "source": [
    "When displaying `pandas` tables in the notebook,\n",
    "only display up to 20 rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14d46dfe-b02c-441e-92fb-f3066daff6d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "pandas.set_option('display.max_rows', 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad0d6bef-485d-43bd-9ea2-cc79d5579cca",
   "metadata": {},
   "source": [
    "The Rubin cutout service allows the user to save cutouts as fits files locally on the Rubin Science Platform (RSP).\n",
    "Create a new temporary folder in the home directory called \"temp\" in which to save these files.\n",
    "At the end of the notebook, the last cell will clear the contents and remove the temp folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "155fd661-08f8-4923-8376-5d6fe5334f4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tempdir = 'dp02_13a_temp'\n",
    "temppath = os.getenv(\"HOME\") + '/' + tempdir\n",
    "if not os.path.exists(temppath):\n",
    "    os.makedirs(temppath)\n",
    "    print('Created: ', temppath)\n",
    "else:\n",
    "    print('Already exists: ', temppath)\n",
    "del temppath"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff95dea9-1223-4573-adce-bf85285dd2a0",
   "metadata": {},
   "source": [
    "The following cells define a number of functions to be used throughout the tutorial. The first function is to enable easy plotting of files from the cutout tool (`plotImage`). The next function, `make_image_cutout`, is a wrapper function that will perform the operations needed to call the cutout tool and create image cutouts stored locally. This procedure will first be demonstrated in Section 2. The steps include: 1) define the location on the sky. 2) Query the TAP service for the specifications of the dataId. 3) Retrieve the datalink URL associated with the data. 4) Create a cutout instance from the query result and the Datalink Resource."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f464394-534d-4164-93d4-eaa27f46bae5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotImage(exposure: ExposureF):\n",
    "    \"\"\"Plot and image using lsst.awf.image package\n",
    "\n",
    "   Parameters\n",
    "    ----------\n",
    "    exposure : `ExposureF`\n",
    "        the image to plot from file in LSST awf image exposure class format\n",
    "\n",
    "   Returns\n",
    "    -------\n",
    "    image for notebook display\n",
    "    \"\"\"\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    display = afwDisplay.Display(frame=fig)\n",
    "    display.scale('asinh', 'zscale')\n",
    "    display.mtv(exposure.image)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16b8a3f9-c6fb-404d-a07c-c852e6afc645",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def make_image_cutout(tap_service, ra, dec, dataId, cutout_size=0.01,\n",
    "                      imtype=None, filename=None):\n",
    "    \"\"\"Wrapper function to generate a cutout using the cutout tool\n",
    "\n",
    "   Parameters\n",
    "    ----------\n",
    "    tap_service : an instance of the TAP service\n",
    "    ra, dec : 'float'\n",
    "        the ra and dec of the cutout center\n",
    "    dataId : 'dict'\n",
    "        the dataId of the image to make a cutout from. The format\n",
    "        must correspond to that provided for parameter 'imtype'\n",
    "    cutout_size : 'float', optional\n",
    "        edge length in degrees of the cutout\n",
    "    imtype : 'string', optional\n",
    "        string containing the type of LSST image to generate\n",
    "        a cutout of (e.g. deepCoadd, calexp). If imtype=None,\n",
    "        the function will assume a deepCoadd.\n",
    "    filename : 'string', optional\n",
    "        filename of the resulting cutout (which has fits format)\n",
    "\n",
    "   Returns\n",
    "    -------\n",
    "    sodaCutout : 'string'\n",
    "        filename of the cutout in fits format (including\n",
    "        full directory path; for now by default it is saved\n",
    "        in /home/dp02_13a_temp/\n",
    "    \"\"\"\n",
    "\n",
    "    spherePoint = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)\n",
    "\n",
    "    if imtype == 'calexp':\n",
    "\n",
    "        query = \"SELECT access_format, access_url, dataproduct_subtype, \" + \\\n",
    "            \"lsst_visit, lsst_detector, lsst_band \" + \\\n",
    "            \"FROM ivoa.ObsCore WHERE dataproduct_type = 'image' \" + \\\n",
    "            \"AND obs_collection = 'LSST.DP02' \" + \\\n",
    "            \"AND dataproduct_subtype = 'lsst.calexp' \" + \\\n",
    "            \"AND lsst_visit = \" + str(dataId[\"visit\"]) + \" \" + \\\n",
    "            \"AND lsst_detector = \" + str(dataId[\"detector\"])\n",
    "        results = tap_service.search(query)\n",
    "\n",
    "    else:\n",
    "        # Find the tract and patch that contain this point\n",
    "        tract = dataId[\"tract\"]\n",
    "        patch = dataId[\"patch\"]\n",
    "\n",
    "        # add optional default band if it is not contained in the dataId\n",
    "        band = dataId[\"band\"]\n",
    "\n",
    "        query = \"SELECT access_format, access_url, dataproduct_subtype, \" + \\\n",
    "            \"lsst_patch, lsst_tract, lsst_band \" + \\\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(tract) + \" \" + \\\n",
    "            \"AND lsst_patch = \" + str(patch) + \" \" + \\\n",
    "            \"AND lsst_band = \" + \"'\" + str(band) + \"' \"\n",
    "        results = tap_service.search(query)\n",
    "\n",
    "    # Get datalink\n",
    "    dataLinkUrl = results[0].getdataurl()\n",
    "    auth_session = service._session\n",
    "    dl = DatalinkResults.from_result_url(dataLinkUrl,\n",
    "                                         session=auth_session)\n",
    "\n",
    "    # from_resource: creates a instance from\n",
    "    # a number of records and a Datalink Resource.\n",
    "    sq = SodaQuery.from_resource(dl,\n",
    "                                 dl.get_adhocservice_by_id(\"cutout-sync\"),\n",
    "                                 session=auth_session)\n",
    "\n",
    "    sq.circle = (spherePoint.getRa().asDegrees() * u.deg,\n",
    "                 spherePoint.getDec().asDegrees() * u.deg,\n",
    "                 cutout_size * u.deg)\n",
    "\n",
    "    tempdir = 'dp02_13a_temp/'\n",
    "\n",
    "    if filename:\n",
    "        sodaCutout = os.path.join(os.getenv('HOME'), tempdir+filename)\n",
    "    else:\n",
    "        sodaCutout = os.path.join(os.getenv('HOME'), tempdir+'cutout.fits')\n",
    "\n",
    "    with open(sodaCutout, 'bw') as f:\n",
    "        f.write(sq.execute_stream().read())\n",
    "\n",
    "    return sodaCutout"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11b72063-a314-409b-8fe2-28d909dfcf6b",
   "metadata": {},
   "source": [
    "Finally, this third cell below defines a set of functions that we will use to work with cutouts of the same source to demonstrate how to visualize variability between calexp files, for diaObjects. Since calexp files may be imaged at any orientation on the sky, in order to align the images the user will 1) rotate and scale the image to a common grid (`warp_image`) and 2) determine the bounding box of the cutouts (`get_minmax_xy`) and 3) create a gif of those warped images to visualize the variability of the source (`make_gif`). These will be demonstrated in Section 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f27addf2-cb60-4456-a6f4-8b8d0bc1b74e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def warp_img(ref_img, img_to_warp, ref_wcs, wcs_to_warp):\n",
    "    \"\"\"Warp and rotate an image onto the coordinate system of another image\n",
    "\n",
    "   Parameters\n",
    "    ----------\n",
    "    ref_img: 'ExposureF'\n",
    "        is the reference image for the re-projection\n",
    "    img_to_warp: 'ExposureF'\n",
    "        the image to rotate and warp onto the reference image's wcs\n",
    "    ref_wcs: 'wcs object'\n",
    "        the wcs of the reference image (i.e. ref_img.getWcs() )\n",
    "    wcs_to_warp: 'wcs object'\n",
    "        the wcs of the image to warp (i.e. img_to_warp.getWcs() )\n",
    "   Returns\n",
    "    -------\n",
    "    warpedExp: 'ExposureF'\n",
    "        a reprojected, rotated image that is aligned and matched to ref_image\n",
    "     \"\"\"\n",
    "    config = RegisterConfig()\n",
    "    task = RegisterTask(name=\"register\", config=config)\n",
    "    warpedExp = task.warpExposure(img_to_warp, wcs_to_warp, ref_wcs,\n",
    "                                  ref_img.getBBox())\n",
    "\n",
    "    return warpedExp\n",
    "\n",
    "\n",
    "def get_minmax_xy(img, cutout_size):\n",
    "    \"\"\"Get the pixel dimensions of an image\n",
    "\n",
    "   Parameters\n",
    "    ----------\n",
    "    img: 'ExposureF'\n",
    "        is the input image to return the pixel coordinates of\n",
    "    cutout_size: 'int'\n",
    "        is the edge size of the image in pixels\n",
    "   Returns\n",
    "    -------\n",
    "    minx, maxx, miny, maxy: 'int'\n",
    "        the min and max x and y pixel coordinates for the input image\n",
    "     \"\"\"\n",
    "\n",
    "    cutout_size = int(cutout_size)\n",
    "\n",
    "    height = img.height\n",
    "    width = img.width\n",
    "\n",
    "    ceny = (height - 1) / 2\n",
    "    cenx = (width - 1) / 2\n",
    "\n",
    "    minx = int(cenx - ((cutout_size - 1) / 2))\n",
    "    maxx = int(cenx + ((cutout_size - 1) / 2))\n",
    "    miny = int(ceny - ((cutout_size - 1) / 2))\n",
    "    maxy = int(ceny + ((cutout_size - 1) / 2))\n",
    "\n",
    "    return {'minx': minx, 'maxx': maxx, 'miny': miny, 'maxy': maxy}\n",
    "\n",
    "\n",
    "def make_gif(folder):\n",
    "    \"\"\"Generate a GIF for all *.png images contained in a folder\n",
    "\n",
    "   Parameters\n",
    "    ----------\n",
    "    # folder: 'string'\n",
    "        string containing the path to the folder\n",
    "        default filename is animation.gif\n",
    "\n",
    "   Returns\n",
    "    -------\n",
    "     \"\"\"\n",
    "    frames = [Image.open(img) for img in sorted(glob.glob(f\"{folder}/*.png\"))]\n",
    "    frame_one = frames[0]\n",
    "    frame_one.save(folder+\"/animation.gif\", format=\"GIF\",\n",
    "                   append_images=frames, save_all=True, duration=500, loop=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec51ac0b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. A step by step demonstration of how to use the Rubin Image Cutout Service\n",
    "\n",
    "This section will demonstrate a simple example: how to create a cutout for a single part of sky from a deepCoadd. \n",
    "\n",
    "\n",
    "### 2.1 Initiate the TAP service, and define sky coordinates for the image cutout\n",
    "\n",
    "The TAP service is used to query the ivoa.Obscore table for the datalink (a web URL identifying where the data is hosted)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f632c1bc-44dd-4dff-9f29-3e63bf89428d",
   "metadata": {},
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47443c3-9455-4707-903f-59ba991cb014",
   "metadata": {},
   "source": [
    "First, define a point on the sky as the center of the image cutout. This example uses the galaxy cluster from DP0.2 Notebook Tutorial 03a. Once the RA and Dec are defined, create a SpherePoint object to define the location on the sky. Use these in a TAP query to identify the overlapping deepCoadd image.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc11b12-7c0b-4ab8-88b0-06787ffd9fe8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ra = 55.7467\n",
    "dec = -32.2862\n",
    "coord = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')\n",
    "radius = .1 * u.deg\n",
    "\n",
    "spherePoint = lsst.geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "551c859a-730a-4f4e-811b-6c2a13c4c30d",
   "metadata": {},
   "source": [
    "The next cell shows the TAP query for the metadata that is associated with the image's remote location in the LSST data archive. The DP0.2 has a schema (table collection) called \"ivoa\", which contains a table called ivoa.ObsCore. The IVOA-defined obscore table contains generic metadata for the DP0.2 datasets. The table is accessible via ADQL queries via a TAP endpoint. The mechanism for locating the images LSST is to use the TAP service to query the ObsCore schema. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80d3f784-e547-4834-b5a4-37381bb15481",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT lsst_patch, lsst_tract, s_region, \" + \\\n",
    "        \"access_format, access_url, dataproduct_subtype \" + \\\n",
    "        \"FROM ivoa.ObsCore \" + \\\n",
    "        \"WHERE dataproduct_type = 'image' \" + \\\n",
    "        \"AND obs_collection = 'LSST.DP02' \" + \\\n",
    "        \"AND dataproduct_subtype = 'lsst.deepCoadd_calexp' \" + \\\n",
    "        \"AND lsst_band = 'i' \" + \\\n",
    "        \"AND CONTAINS(POINT('ICRS', \" + str(coord.ra.value) + \\\n",
    "        \", \" + str(coord.dec.value) + \"), \" + \\\n",
    "        \"s_region) = 1\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "844332b9-17a3-420b-b410-034bda4be211",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a97c1a71-58c3-45f6-98dc-9933299b5444",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcf13212-6207-4336-9b3f-e0f7aaaaa6a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "317b66e1-8a2b-4268-803d-08c1cb168350",
   "metadata": {},
   "source": [
    "In the above table, the `access_url` contains the web URL datalink for the requested deepCoadd. This datalink will be needed to generate the image cutout.\n",
    "\n",
    "The identified Patch and Tract will be used to define the unique dataId for that location, and once a filter (band) is included, this defines a unique deepCoadd in the LSST image database. Below, construct the dataId."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f075025-e609-4c1c-a965-eda296fcd549",
   "metadata": {},
   "outputs": [],
   "source": [
    "tract = results['lsst_tract'][0]\n",
    "patch = results['lsst_patch'][0]\n",
    "\n",
    "dataId = {'band': 'i', 'tract': tract,\n",
    "          'patch': patch}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "996dd37e-7c19-4796-8b19-5362d347bc29",
   "metadata": {},
   "source": [
    "### 2.2 Generating an image cutout \n",
    "\n",
    "First, extract the datalink (`dataLinkUrl`) from the query result using the `getdataurl` method. Next, extract the session authentication, which is needed for reuse of the datalink. Finally, create a `DatalinkResults` object to be able to access this URL, which will be stored as `dl_results` and available for approximately 15 minutes, in a format that can be used by the IVOA tools below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb4b029f-f4b5-4f51-8483-96608c2711ae",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dataLinkUrl = results[0].getdataurl()\n",
    "\n",
    "auth_session = service._session\n",
    "dl_results = DatalinkResults.from_result_url(dataLinkUrl,\n",
    "                                             session=auth_session)\n",
    "\n",
    "f\"Datalink status: {dl_results.status}. Datalink service url: {dataLinkUrl}\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb03b3a3-d067-4e36-bc9b-754d966c1477",
   "metadata": {},
   "source": [
    "Lastly, call the Rubin Image Cutout Service. It is done by feeding the data link created above (called `dl_results`) to following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c667776-5f21-4879-b743-89b89b78337d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = SodaQuery.from_resource(dl_results,\n",
    "                             dl_results.get_adhocservice_by_id(\"cutout-sync\"),\n",
    "                             session=auth_session)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac93f7d9-dca1-43dc-b01c-1d1b4b5d145b",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2023-06-01T21:40:00.560310Z",
     "iopub.status.idle": "2023-06-01T21:40:00.560562Z",
     "shell.execute_reply": "2023-06-01T21:40:00.560451Z",
     "shell.execute_reply.started": "2023-06-01T21:40:00.560440Z"
    }
   },
   "source": [
    "The variable `sq` now holds the result of the SODA query using the data link (which currently still points the full LSST patch and tract deepCoadd image, at its remote location in the database). The cell below will now demonstrate how to extract a cutout from `sq`.  \n",
    "\n",
    "### 2.2.1 Defining a square cutout using a subtended circle\n",
    "\n",
    "Only two shape definitions are supported: a circle function, and a polygon function can be used to define the cutout dimensions. These shape definitions do not produce circle or polygon cutouts, but rather are methods for defining the edges of cutouts with 4 sides. In the case of circle, the resulting cutout is always a square, with edge size that is the same as the circle diameter. In the case of a polygon, either a square or a rectangular cutout will result, depending on whether the length and width edge dimensions are different values. Only cutouts with 4 corners and 90 degree angles are supported. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23ca3f11-f17e-4e5c-8469-a0d25554234d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sphereRadius = 0.05 * u.deg\n",
    "sq.circle = (spherePoint.getRa().asDegrees() * u.deg,\n",
    "             spherePoint.getDec().asDegrees() * u.deg,\n",
    "             sphereRadius)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e3e8f5-bcde-4b87-85ab-1acdd00ebce6",
   "metadata": {},
   "source": [
    "### 2.2.2 Saving and plotting the cutout\n",
    "\n",
    "This section will demonstrate two ways to handle the output of the cutout service. While the image cutout procedure itself is performed remotely, bytes are returned, and the user can decide how to receive it. The next cell demonstrates how to create a fits file, stored in the returned Google cloud storage service (GCS) bucket. In the following example, the cutout is written to disk as a fits file (which is done using the write() function), before it is accessed by the notebook. (Future planned options for the Rubin cutout service, including the potential to retrieve other image formats such as jpeg, are listed at the <a href=\"https://dmtn-208.lsst.io/\">Rubin Science Platform image cutout implementation strategy document</a>).\n",
    "\n",
    "All cutouts saved to disk will be stored in the temporary folder in this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12183fdf-7d11-4268-b7ac-025a9e904230",
   "metadata": {},
   "outputs": [],
   "source": [
    "sodaCutout = os.path.join(os.getenv('HOME'), tempdir + '/cutout-circle.fits')\n",
    "with open(sodaCutout, 'bw') as f:\n",
    "    f.write(sq.execute_stream().read())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13e154ac-ea08-49bf-8799-192de500e71b",
   "metadata": {},
   "source": [
    "Finally, call the `plotImage` function (defined in Section 1), which will quickly display the cutout that was written to disk in in the temporary folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d88e832b-508c-4d29-8748-320c2d1fded9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotImage(ExposureF(sodaCutout))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "408e0448-59f8-4db9-b4f6-ed44d11496ee",
   "metadata": {},
   "source": [
    "> Figure 1: The `deepCoadd` cutout image, displayed in grayscale with a scale bar at right, and image axes in pixels. The field is centered on a rich galaxy cluster.\n",
    "\n",
    "Option: the cutout can be stored directly in memory without first writing to disk. This method is demonstrated in the below cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea228c0-ae77-4e65-aceb-2bec77d758d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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",
    "# plotImage(exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99e35e80-15ce-40c6-aa8b-9079e3bc0200",
   "metadata": {},
   "source": [
    "### 2.2.3 Using polygon to define the image cutout shape instead of a circle\n",
    "\n",
    "It is also possible to define the cutout geometry using a polygon, which enables the cutout to be rectangular, but not necessarily be square. For this, use `polygon`, which takes as input the four corners in celestial coordinates. A minimum of 3 vertices are required (the line from the last vertex back to the first is implicit) Vertices must be ordered in the counter-clockwise direction. For example: a polygon is defined as a set of 4 (x,y) coordinates from (12,34) to (14,34) to (14,36) to (12,36) and (implicitly) back to (12,34) as:\n",
    "\n",
    "POLYGON=12 34 14 34 14 36 12 36\n",
    "\n",
    "Since the center of the galaxy cluster is already defined in RA and Dec in the cells above (`spherePoint`), this example will define each x,y set as RA+/-sphereRadius and Dec+/-sphereRadius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef262797-3af1-48fb-955a-7a3d2ebe4c55",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq2 = SodaQuery.from_resource(dl_results,\n",
    "                              dl_results.get_adhocservice_by_id(\"cutout-sync\"),\n",
    "                              session=auth_session)\n",
    "sphereRadius1 = 0.03 * u.deg\n",
    "sphereRadius2 = 0.01 * u.deg\n",
    "\n",
    "sq2.polygon = (spherePoint.getRa().asDegrees() * u.deg - sphereRadius1,\n",
    "               spherePoint.getDec().asDegrees() * u.deg - sphereRadius2,\n",
    "               spherePoint.getRa().asDegrees() * u.deg - sphereRadius1,\n",
    "               spherePoint.getDec().asDegrees() * u.deg + sphereRadius2,\n",
    "               spherePoint.getRa().asDegrees() * u.deg + sphereRadius1,\n",
    "               spherePoint.getDec().asDegrees() * u.deg + sphereRadius2)\n",
    "\n",
    "sodaCutout2 = os.path.join(os.getenv('HOME'), tempdir + '/cutout-polygon.fits')\n",
    "with open(sodaCutout2, 'bw') as f:\n",
    "    f.write(sq2.execute_stream().read())\n",
    "\n",
    "plotImage(ExposureF(sodaCutout2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6970143e-90b9-406c-80ce-e6c6ba923c1e",
   "metadata": {},
   "source": [
    "> Figure 2: A rectangular polygon cutout of the image displayed in Figure 1.\n",
    "\n",
    "There is an important difference to note between the circle and polygon shape definitions. The angular distance on the sky that defines the circular cutout size already accounts for the difference in angular distance in the RA direction is smaller by a factor of cos(declination), where declination is in units radians. The difference increases with higher declination. However, the polygon definition does not automatically account for this cosine factor. Thus, circle and polygon cutout definitions using the same cutout edge length will not match size in the RA direction (for deepCoadds, the x-axis). The cell below demonstrates how to make this correction to the polygon cutout definition to create symmetric cutouts with polygon. Here `sphereRadius` comes from the circle definition above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3407da92-6494-4ef0-9343-e5f4db20a87d",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = Angle(spherePoint.getDec().asDegrees(), u.deg)\n",
    "cosd = np.cos(a.radian)\n",
    "\n",
    "sq2.polygon = (spherePoint.getRa().asDegrees() * u.deg - sphereRadius/cosd,\n",
    "               spherePoint.getDec().asDegrees() * u.deg - sphereRadius,\n",
    "               spherePoint.getRa().asDegrees() * u.deg - sphereRadius/cosd,\n",
    "               spherePoint.getDec().asDegrees() * u.deg + sphereRadius,\n",
    "               spherePoint.getRa().asDegrees() * u.deg + sphereRadius/cosd,\n",
    "               spherePoint.getDec().asDegrees() * u.deg + sphereRadius,\n",
    "               spherePoint.getRa().asDegrees() * u.deg + sphereRadius/cosd,\n",
    "               spherePoint.getDec().asDegrees() * u.deg - sphereRadius)\n",
    "\n",
    "sodaCutout2 = os.path.join(os.getenv('HOME'), tempdir + '/cutout_polygon.fits')\n",
    "with open(sodaCutout2, 'bw') as f:\n",
    "    f.write(sq2.execute_stream().read())\n",
    "\n",
    "plotImage(ExposureF(sodaCutout2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3347e8f2-fb13-44fb-92d9-e6d6c64ad8be",
   "metadata": {},
   "source": [
    "> Figure 3: The same cutout as displayed in Figure 1, but using the `polygon` functionality and accounting for the $cos({\\rm dec})$ term."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c53ba0bc-91c2-4576-8014-4e008929d50c",
   "metadata": {},
   "source": [
    "### 2.3 Test out the cutout wrapper function \"make_image_cutout\" \n",
    "\n",
    "All of the above steps from section 2.2 have been compiled in a wrapper function called `make_image_cutout` that is defined in Section 1. In the rest of this notebook, cutouts will be generated using this function for simplicity. As defined, it assumes the circular cutout definition demonstrated above. Thus, the function requires as input the TAP service, the center ra/dec of the cutout, the `dataId` and imtype (Section 4 will demonstrate how to do this for calexp, not just deepCoadd) and the size of the cutout (i.e. what was defined as `sphereRadius` above). The next cell demonstrates how to run all the steps by calling `make_image_cutout` and plotting the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd1a68b4-5524-4016-a743-88b34e1d661a",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "imtype = 'deepCoadd'\n",
    "test = make_image_cutout(service, ra, dec, dataId=dataId,\n",
    "                         imtype=imtype, cutout_size=0.005)\n",
    "plotImage(ExposureF(test))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cbfc19f-3fd9-4c5e-a3a9-e9ed9f941589",
   "metadata": {},
   "source": [
    "> Figure 4: A small image cutout, created using the `make_image_cutout` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9662577f-594f-4442-9c72-654bc8aef4b7",
   "metadata": {},
   "source": [
    "## 3. Science application: validating a galaxy sample based on color-color selection \n",
    "\n",
    "This section will demonstrate a simple example of a science use-case for the Rubin Image Cutout Service. \n",
    "\n",
    "Lyman-break color selections are a common way to identify galaxies at high redshifts. The selections make use of the absorption by intergalactic hydrogen in the foreground of all galaxy light blueward of the Lyman limit (912A) or at very high redshifts, Lyman-alpha (1216A). <a href=\"http://candels-collaboration.blogspot.com/2012/08/how-to-find-distant-galaxies.html\">More details about why this works can be found at this blog-post.</a> DP0.2 includes mock galaxies up to redshift (z) of 3, where the Lyman-break falls between the LSST u and g bands. Such galaxies can be identified by their very red u-g colors, and since intergalactic hydrogen absorbs all the light blueward of the Lyman break, true z~3 galaxies should be undetected in the LSST u-band. The cell below will perform a search for redshift of 3 galaxies as u-band dropouts, and will create cutouts in each of the LSST filters in order to visually inspect the Lyman-break galaxy (LBG) candidates to confirm the real candidates show no u-band flux. \n",
    "\n",
    "### 3.1 The u-band dropout LBG color selection\n",
    "To search for $z\\sim3$ LBGs in DP0.2, first define a u-band color selection based on the existing literature that uses filter sets similar to that of LSST. Below outlines the set of colors used in Steidel et al. 2003 (see also Boutsia et al 2014, among other references. Many selections exist, where changes to colors would alter the redshift selection function, or may differ due to differences in the broad band filter shapes of other instruments).\n",
    "\n",
    "U − G >= (G − R) + 1.0\n",
    "\n",
    "G − R <= 1.2\n",
    "\n",
    "\n",
    "where capital letters U, G, and R indicate the total AB magnitudes in the LSST filters u, g, r. Since DP0.2 is only 5-year depth mock data in the LSST filters, our search will also include a g-band magnitude that is relatively bright (G < 24.5 ABmag) in order to obtain high-quality candidates. \n",
    "\n",
    "High redshift galaxies are typically small in size. For this exercise, use aperture photometry with a relatively small aperture diameter (9 pixels) in order to obtain high S/N measurements and exclude neighbors.\n",
    "\n",
    "#### 3.1.1 Ammending a typical search because of the rarity of z~3 galaxies in the mock data.\n",
    "The mock galaxies that go into DP0.2 are created using the LSST Catalog Simulator (<a href=https://www.lsst.org/scientists/simulations/catsimcatSim>catSim</a>), and are based on an empirical model of the evolution of galaxies outlined in <a href=\"https://arxiv.org/abs/1907.06530\">Korytov et al. 2019</a>. The input model is based on realistic number densities, fluxes, and redshift distributions across cosmic time, and utilizes the UniverseMachine methodology to assign empirical properties (see <a href=\"https://arxiv.org/abs/1806.07893\">Behroozi et al. 2019</a>). \n",
    "\n",
    "Typically at any redshift, bright galaxies are rarer than faint galaxies, which are much more common. And typically in any patch of the sky, nearby galaxies will be much more common above the detection limit of a survey than distant $z\\sim3$ galaxies. This means that $z\\sim3$ galaxies that are bright enough to meet the criteria we want (G$ < 24.5$) will require a large search area to identify a statistical sample. So, in this example, we will first identify a parent sample of bright (G < 24.5) galaxies over a very large area (4 degrees) to ensure the parent sample contains enough rare distant galaxies.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5c745ad-3db2-47cc-b50b-5d0beae29a40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "center_coords = \"62, -37\"\n",
    "ralbg = 62.\n",
    "declbg = -37.\n",
    "radius = \"4.0\"\n",
    "max_rec = 500000\n",
    "\n",
    "query = \"SELECT TOP \" + str(max_rec) + \" \" + \\\n",
    "        \"objectId, coord_ra, coord_dec, detect_isPrimary, patch, tract, \" + \\\n",
    "        \"u_ap09Flux, u_ap09FluxErr, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(u_ap09Flux) as umag, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(r_ap09Flux) as rmag, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(g_ap09Flux) as gmag, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(u_ap09FluxErr) as umagErr \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + center_coords + \", \" + radius + \")) = 1 \" + \\\n",
    "        \"AND detect_isPrimary = 1 AND detect_fromBlend = 0 \" + \\\n",
    "        \"AND g_ap09Flux/g_ap09FluxErr > 10 AND \" + \\\n",
    "        \"r_ap09Flux/r_ap09FluxErr > 10\" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(g_ap09Flux) < 24.5 \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19bb33f2-5257-46fe-a333-eb5dcb7b20eb",
   "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,
   "id": "c6be498e-871b-4d45-bfad-bd89ca607106",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c3e8bae-e010-4b2e-b3fc-e688ebf42c5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77189c9e-85b8-441f-9d93-ba1277bb598d",
   "metadata": {},
   "source": [
    "The next cells create a plot of the U-G vs G-R colors of the parent sample, and identify the high-redshift galaxy candidates based on the color selection discussed above. For sources with low S/N fluxes in u-band, we cannot robustly recover the U-G color, and instead we get a lower limit on how red it is (owing to the upper limit to the u-band flux). Thus, the U-G color can be overestimated if the u-band flux is not significant. Thus, sources with low S/N u-band flux will be reset to an error floor that corresponds the noise level of the u-band imaging so that it's color is not overestimated. To do this, define a new parameter called `umag` to store the u-band magnitude for assessing colors. In this case, if sources have signal to noise less than 1 in U_band, umag will instead hold the `umagErr` from the object table)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5627dc74-9e03-4dd7-b6a3-fa328b7ee9cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "umag = np.zeros(len(results['umag']), dtype='float')\n",
    "S2N = np.asarray(results['u_ap09Flux']/results['u_ap09FluxErr'], dtype='float')\n",
    "\n",
    "for i in range(len(umag)):\n",
    "    if S2N[i] < 1:\n",
    "        umag[i] = results['umagErr'][i]\n",
    "    else:\n",
    "        umag[i] = results['umag'][i]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "937a41f4-d45c-45a7-bb6f-ef24dd506c2e",
   "metadata": {},
   "source": [
    "Further, this simple example uses a restricted color selection that is more narrow and blue in G-R than typical selections (G-R$<$0.2; this may exclude real z=3 galaxies that are redder in G-R, for example galaxies that have older stars or are dusty; but the restriction is helpful for this simple example because it lowers the number of interloper sources). The LBG candidates (indexed from the array as `whlbg`) identified using the color selection (black dashed lines) are in orange, with the full query of all objects in blue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51f4bb38-91b5-4caf-b92c-4ba1319a307e",
   "metadata": {},
   "outputs": [],
   "source": [
    "GRlim = 0.2\n",
    "GRfloor = 1.2\n",
    "\n",
    "whlbg = np.where(((umag-results['gmag'])\n",
    "                  >= ((results['gmag'] - results['rmag']) + GRfloor))\n",
    "                 & ((results['gmag'] - results['rmag']) < GRlim)\n",
    "                 & ((results['gmag'] - results['rmag']) > 0))[0]\n",
    "\n",
    "GmR = np.arange(0, GRlim + 0.01, 0.1)\n",
    "plt.plot(GmR, GmR + GRfloor, '--', color='k')\n",
    "plt.plot([0, 0], [GRfloor, 3], '--', color='k')\n",
    "plt.plot([GRlim, GRlim], [GRlim + GRfloor, 3], '--', color='k')\n",
    "\n",
    "GminusR = results['gmag']-results['rmag']\n",
    "UminusG = umag-results['gmag']\n",
    "\n",
    "plt.plot(GminusR, UminusG, '.', label='all objects', alpha=.1)\n",
    "plt.plot(GminusR[whlbg], UminusG[whlbg], '.',\n",
    "         label='color-selected LBG candidates')\n",
    "\n",
    "plt.xlabel('G - R')\n",
    "plt.ylabel('U - G')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f93f4a8a-d624-4d68-8591-0ebe4d0c5a8f",
   "metadata": {},
   "source": [
    "> Figure 5: The $g-r$ versus $u-g$ color of all objects (blue) and color-selected LBG candidates (orange), with the box defining LBG candidates marked with black dashed lines.\n",
    "\n",
    "Some fraction of the orange LBG candidates might be interlopers, entering the selection box for various reasons (e.g. photometric scatter). The next two cells below will plot image cutouts of each LBG candidate for visual inspection. Since this is a simulated dataset, the first cell performs `query3` to the TAP service to pull the intrinsic redshift `ts_redshift` for each object from the truth table. This will be used to compare with the expected redshift $z\\sim3$ of selected objects in the image cutout captions (see DP0.2 Tutorial Notebook 08 for more information on truth tables).\n",
    "\n",
    "To limit the compute time (dominated by creating the cutouts), select just the first 5 galaxy candidates among the sample. This number can be altered using `n_plot`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45b2ee57-7d00-4b3a-a485-2572a011be96",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_plot = 5\n",
    "\n",
    "id = str(tuple(list(results['objectId'][whlbg][0:n_plot])))\n",
    "\n",
    "query3 = \"SELECT mt.id_truth_type AS mt_id_truth_type, \" + \\\n",
    "         \"obj.objectId AS obj_objectId,  \" + \\\n",
    "         \"ts.redshift AS ts_redshift \" + \\\n",
    "         \"FROM dp02_dc2_catalogs.MatchesTruth AS mt \" + \\\n",
    "         \"JOIN dp02_dc2_catalogs.TruthSummary AS ts \" + \\\n",
    "         \"ON mt.id_truth_type=ts.id_truth_type \" + \\\n",
    "         \"JOIN dp02_dc2_catalogs.Object AS obj \" + \\\n",
    "         \"ON mt.match_objectId=obj.objectId \" + \\\n",
    "         \"WHERE obj.objectId IN \"+id+\" \" + \\\n",
    "         \"AND ts.truth_type=1 \" + \\\n",
    "         \"AND obj.detect_isPrimary=1 \"\n",
    "\n",
    "results3 = service.search(query3)\n",
    "results3.to_table()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77bda4aa-398c-47e3-957c-0d77c7b35f4c",
   "metadata": {},
   "source": [
    "The cell below creates the cutouts for each of the objects in 5 different bands. By inspecting the cutouts and intrinsic redshift, there are a small number of interlopers (at redshifts $z<2.5$ or so) which visibly show some u-band flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "799441f5-cb29-42a9-94cb-9c345ff1818f",
   "metadata": {},
   "outputs": [],
   "source": [
    "filt = [\"u\", \"g\", \"r\", \"i\", \"z\"]\n",
    "\n",
    "for i in range(n_plot):\n",
    "    fig, axs = plt.subplots(1, len(filt), figsize=(25, 25))\n",
    "\n",
    "    for j, ax in enumerate(fig.axes):\n",
    "        dataId_deep = {'patch': results['patch'][whlbg][i],\n",
    "                       'tract': results['tract'][whlbg][i],\n",
    "                       'band': filt[j]}\n",
    "        cutout = make_image_cutout(service, results['coord_ra'][whlbg][i],\n",
    "                                   results['coord_dec'][whlbg][i],\n",
    "                                   cutout_size=0.001, imtype='deepCoadd',\n",
    "                                   dataId=dataId_deep,\n",
    "                                   filename='cutout_' + filt[j] + '_'\n",
    "                                   + str(i) + '.fits')\n",
    "\n",
    "        hdulist = fits.open(cutout)\n",
    "        img_hdr = hdulist[1].header\n",
    "        img_wcs = WCS(img_hdr)\n",
    "        img_data = fits.getdata(cutout)\n",
    "        ax.imshow(img_data, cmap='gray',\n",
    "                  vmin=0, vmax=0.7, norm='linear',\n",
    "                  origin='lower')\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        ax.text(.1, .9, filt[j]+'-band z=' + str(results3['ts_redshift'][i]),\n",
    "                color='white')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5fa95dd-3323-469f-9c7a-8fff421bdfbf",
   "metadata": {},
   "source": [
    "> Figure 6: Cutout images in $ugriz$ bands (left to right) for 5 of the LBG candidates (top to bottom)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5111d5-82a2-43ec-b4a8-daafdda10e72",
   "metadata": {},
   "source": [
    "## 4. Creating calexp cutouts to inspect changing flux for the diaObject\n",
    "\n",
    "Another useful application of the cutout tool is to visually inspect a large number of images of a source identified with difference image analysis (DIA) in order to confirm its variability. \n",
    "\n",
    "This example is different in that it requires generating cutouts of calexp files around the time the transient is identified, rather than the deepCoadds as used above. \n",
    "\n",
    "### 4.1 Query for the source information needed to plot images before and after the transient event\n",
    "\n",
    "This example will use a known diaObject with ID 1250953961339360185 (one identified and studied in DP0.2 tutorial notebook 07a). The cell below will retrieve the information needed to generate a set of images in order of time, including and after the transient event. The TAP query will identify 8 diaSource IDs measured for this diaObject at the ra/dec, and will also retrieve the ccdVisitId of each time the diaSource was observed in a calexp image.  \n",
    "\n",
    "The last line will sort the query return by date, so that they appear in the correct chronological order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d39ecde-b94c-4cff-b6bc-3219d50ecebf",
   "metadata": {},
   "outputs": [],
   "source": [
    "diaObjectId = 1250953961339360185\n",
    "query = \"SELECT ra, decl, diaObjectId, diaSourceId, \"\\\n",
    "        \"psFlux, psFluxErr, filterName, midPointTai, \"\\\n",
    "        \"SNR, ccdVisitId, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(psFlux) AS psAbMag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaSource \"\\\n",
    "        \"WHERE diaObjectId = \" + str(diaObjectId)\\\n",
    "        + \" AND filterName = 'i' \"\\\n",
    "        \"ORDER BY midPointTai\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2847df68-19eb-4113-85b3-c96a2e1d1a61",
   "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,
   "id": "40b26604-e91e-4a2b-acb6-5a9b06bbd3e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "diaSrc = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b1b63af-fead-4e3f-8d75-af70ed2bae4b",
   "metadata": {},
   "source": [
    "The next cell will retrieve calexp cutouts for each of the dia sources in the query, and plot them. Vmin and vmax set a common scaling so that the image scale is not skewed by a bright neighbor, and emphasizes the lower brightness sources such as the dia source. Each figure also is labeled by the time since the first difference image ($\\Delta$t in units MJD)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a5bdf81-8438-4cb8-806b-d4e0778da92c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "ra = diaSrc['ra'][0] \n",
    "dec = diaSrc['decl'][0]\n",
    "vmin = -200\n",
    "vmax = 300    \n",
    "i=0\n",
    "for src in diaSrc:\n",
    "    ccdvisitID = src['ccdVisitId']\n",
    "    visit = str(ccdvisitID)[:-3]\n",
    "    detector = str(ccdvisitID)[-3:]\n",
    "    visit = int(visit)\n",
    "    detector = int(detector)\n",
    "\n",
    "    dataId_calexp = {'visit':visit, 'detector':detector}\n",
    "    test = make_image_cutout(service, ra, dec, cutout_size=0.007,\n",
    "                             imtype='calexp', dataId=dataId_calexp,\n",
    "                             filename='cutout_'+str(i)+'.fits')\n",
    "\n",
    "    ref_name = 'cutout_'+str(i)+'.fits'\n",
    "    image_file = os.path.join(os.getenv('HOME'), tempdir + '/' + ref_name)\n",
    "\n",
    "    fig, ax1 = plt.subplots(1,1,figsize=(2, 2))\n",
    "\n",
    "    hdulist = fits.open(image_file)\n",
    "    img_hdr = hdulist[1].header\n",
    "    img_wcs = WCS(img_hdr)\n",
    "    img_data = fits.getdata(image_file)\n",
    "    ax1.imshow(img_data, cmap='gray', \n",
    "               vmin=vmin, vmax=vmax,\n",
    "               origin='lower')\n",
    "    ax1.set_xticks([])\n",
    "    ax1.set_yticks([])\n",
    "    ax1.text(.1, .9,'$\\Delta$t=' + str(round(src['midPointTai']-diaSrc['midPointTai'][0],3)),\n",
    "                color='white',fontsize=12)\n",
    "    plt.show()\n",
    "    i = i+1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "984e239e-6990-4f07-a8dd-47ad00172c14",
   "metadata": {},
   "source": [
    "> Figure 7: The eight individual $i$-band cutout images for the `diaObject`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4197282-3351-4bec-a310-a1f6bc211666",
   "metadata": {},
   "source": [
    "### 4.2 Align the calexp images of the transient\n",
    "\n",
    "In the above plots, it is clear that the calexp image orientations are not aligned with each other, because the calexp images are not necessarily all taken in the same position angle on the sky. The cell below will go through each dia source, identify the calexp file in which it was identified, and create cutouts (as above) except now they will all be aligned with each other. At the same time as creating these new aligned fits files for each cutout, the cell will create png files that will be useful for visualizing the transient behavior (to be used in Section 4.3).\n",
    "\n",
    "The `cutout_size` variable is in units of degrees, converted from `cutout_size_pix` using the LSST Science Camera platescale of 0.2 arcseconds per pixel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88303981-d84d-45f2-b1ad-6a6a761b3141",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "i = 0\n",
    "cutout_size_pix = 131\n",
    "cutout_size = cutout_size_pix * 0.2 / 3600.0\n",
    "\n",
    "for src in diaSrc:\n",
    "    ccdvisitID = src['ccdVisitId']\n",
    "    visit = str(ccdvisitID)[:-3]\n",
    "    detector = str(ccdvisitID)[-3:]\n",
    "    visit = int(visit)\n",
    "    detector = int(detector)\n",
    "    dataId_calexp = {'visit': visit, 'detector': detector}\n",
    "    if i == 0:\n",
    "\n",
    "        # Make the first image the reference image, that all the\n",
    "        # others images will be warped to align with.\n",
    "        ref_name = 'cutout_ref.fits'\n",
    "        ref_cutout = os.path.join(os.getenv('HOME'), tempdir + '/' + ref_name)\n",
    "\n",
    "        make_image_cutout(service, ra, dec, cutout_size=cutout_size,\n",
    "                          imtype='calexp', dataId=dataId_calexp,\n",
    "                          filename=ref_name)\n",
    "\n",
    "        img_ref = ExposureF(ref_cutout)\n",
    "        minmax = get_minmax_xy(img_ref, cutout_size_pix)\n",
    "        im_arr = img_ref.image.array[minmax['minx']: minmax['maxx'],\n",
    "                                     minmax['miny']: minmax['maxy']]\n",
    "\n",
    "        fig, ax = plt.subplots(1, 1, figsize=(2, 2))\n",
    "        plt.imshow(im_arr, origin='lower', cmap='gray', vmin=vmin, vmax=vmax)\n",
    "\n",
    "        ax.set_xticklabels([])\n",
    "        ax.set_yticklabels([])\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "\n",
    "        ax.text(.1, .9, r'$\\Delta$t='\n",
    "                + str(round(src['midPointTai']-diaSrc['midPointTai'][0], 2)),\n",
    "                color='white', fontsize=12)\n",
    "\n",
    "        figname = os.path.join(os.getenv('HOME'),\n",
    "                               tempdir + '/cutout_' + str(i) + '.png')\n",
    "        if os.path.isfile(figname):\n",
    "            os.remove(figname)\n",
    "        plt.savefig(figname)\n",
    "        plt.show()\n",
    "        plt.close()\n",
    "\n",
    "    else:\n",
    "\n",
    "        fitsname = 'cutout_'+str(i)+'.fits'\n",
    "        make_image_cutout(service, ra, dec, cutout_size=cutout_size,\n",
    "                          imtype='calexp', dataId=dataId_calexp,\n",
    "                          filename=fitsname)\n",
    "\n",
    "        img_cutout = os.path.join(os.getenv('HOME'), tempdir + '/'+fitsname)\n",
    "        img = ExposureF(img_cutout)\n",
    "\n",
    "        img_warped = warp_img(img_ref, img, img_ref.getWcs(), img.getWcs())\n",
    "\n",
    "        fig, ax = plt.subplots(1, 1, figsize=(2, 2))\n",
    "\n",
    "        minmax = get_minmax_xy(img_warped, cutout_size_pix)\n",
    "        im_arr = img_warped.image.array[minmax['minx']: minmax['maxx'],\n",
    "                                        minmax['miny']: minmax['maxy']]\n",
    "\n",
    "        plt.imshow(im_arr, origin='lower', cmap='gray', vmin=vmin, vmax=vmax)\n",
    "\n",
    "        ax.set_xticklabels([])\n",
    "        ax.set_yticklabels([])\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        ax.text(.1, .9, r'$\\Delta$t='\n",
    "                + str(round(src['midPointTai']-diaSrc['midPointTai'][0], 2)),\n",
    "                color='white', fontsize=12)\n",
    "\n",
    "        figname = os.path.join(os.getenv('HOME'),\n",
    "                               tempdir + '/cutout_' + str(i) + '.png')\n",
    "        if os.path.isfile(figname):\n",
    "            os.remove(figname)\n",
    "        plt.savefig(figname)\n",
    "        plt.show()\n",
    "        plt.close()\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b921539-cd37-4cd1-95d4-fb1a5ffca9e2",
   "metadata": {},
   "source": [
    "> Figure 8: Similar to Figure 7, but for a smaller cutout size and aligning all of the cutout images.\n",
    "\n",
    "### 4.3 Make a gif of the aligned calexp cutouts for the transient\n",
    "\n",
    "The cell below will take all the png files created above in the temp folder, and make a gif movie that demonstrates the transient behavior of the diaobject. The gif is by default stored in the same directory as the images it is generated from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78a1cb9e-c0cb-46e3-af34-11360b8a6b1a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "images = os.path.join(os.getenv('HOME'), tempdir + '/')\n",
    "make_gif(images)\n",
    "\n",
    "display(dimg(data=open(images+'animation.gif', 'rb').read()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b70402f-11df-403c-8b78-8b8250913fdc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-18T00:10:40.375415Z",
     "iopub.status.busy": "2024-12-18T00:10:40.374694Z",
     "iopub.status.idle": "2024-12-18T00:10:40.379415Z",
     "shell.execute_reply": "2024-12-18T00:10:40.378733Z",
     "shell.execute_reply.started": "2024-12-18T00:10:40.375386Z"
    }
   },
   "source": [
    "> Figure 9: An animated gif centered on a point source (the transient) that brightens then fades."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de4f00c9-04e8-4604-a9a9-dbb1b42004cd",
   "metadata": {},
   "source": [
    "## 5. Remove the temp folder\n",
    "\n",
    "As mentioned in Section 1.2, empty and remove the temporary folder in the home directory.\n",
    "This will remove all files generated by this notebook, and leave the folder empty.\n",
    "\n",
    "To remove the temp folder, open a terminal and type `cd ~` to go to the home directory, then `rmdir dp02_13a_temp` to remove the folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00d1d56b-3dbd-4bb9-8799-771f474c2725",
   "metadata": {},
   "outputs": [],
   "source": [
    "for filepath in glob.glob(os.path.join(os.getenv('HOME'), tempdir + '/')+\"*\"):\n",
    "    os.remove(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0138cbf3-11d9-465e-9b62-7de23ed2e165",
   "metadata": {},
   "source": [
    "## 6. Exercises for the learner\n",
    "\n",
    "1) Extend the time baseline for the cutouts of the DIA object prior to the transient detection (hint, use the ForcedSourceOnDiaObject table)\n",
    "2) Test the failure modes of the cutout tool: what happens if a cutout is requested around coordinates that do not lie inside the linked image?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}