DP02_13a_Image_Cutout_SciDemo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "### <img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<b>Using the image cutout tool with DP0.2</b> <br>\n",
    "Contact author(s): Christina Williams <br>\n",
    "Last verified to run: 2024-4-25 <br>\n",
    "LSST Science Piplines version: Weekly <i>2024_4</i> <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.  **DOESN'T WORK AT CC-IN2P3 AT THE MOMENT, MISSING OBSCORE TABLE\"**"
   ]
  },
  {
   "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-2.lsst.io\">dp0-2.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 operations (e.g. 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> provides a <a href=\"https://www.ivoa.net/documents/SODA/20170517/REC-SODA-1.0.html\">Server-side Operations for Data Access (SODA)</a> for web data access. The LSST architecture is built to fit into the requirements of the Virtual Observatory, enabling the use of Virtual Observatory 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 demonstration of the technical side of this process using DP0.2 data can be found in (Leanne's companion tutorial when published, link here). 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cddc1458",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:09.666188Z",
     "iopub.status.busy": "2024-04-29T14:27:09.665950Z",
     "iopub.status.idle": "2024-04-29T14:27:16.436639Z",
     "shell.execute_reply": "2024-04-29T14:27:16.435809Z",
     "shell.execute_reply.started": "2024-04-29T14:27:09.666166Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Import general python packages\n",
    "import time\n",
    "import numpy as np\n",
    "import uuid\n",
    "import os\n",
    "import glob\n",
    "import math\n",
    "\n",
    "# LSST package imports\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",
    "\n",
    "# Import the Rubin TAP service utilities\n",
    "from lsst.rsp import get_tap_service\n",
    "from lsst.rsp.utils import get_access_token\n",
    "\n",
    "# Packages for warping and aligning calexp cutouts\n",
    "from lsst.pipe.tasks.registerImage import RegisterConfig, RegisterTask\n",
    "\n",
    "# Pandas tools for working with tables\n",
    "import pandas\n",
    "\n",
    "# Plotting functions\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# PyVO packages\n",
    "import pyvo\n",
    "from pyvo.dal.adhoc import DatalinkResults, SodaQuery\n",
    "\n",
    "# Astropy\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\n",
    "\n",
    "from PIL import Image\n",
    "from IPython.display import Image as dimg\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters\n",
    "\n",
    "\n",
    "#### 1.2.1 Set up  plotting defaults for matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "552c4cac-feca-4d51-a92f-87d5779a3300",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:16.438187Z",
     "iopub.status.busy": "2024-04-29T14:27:16.437641Z",
     "iopub.status.idle": "2024-04-29T14:27:16.464081Z",
     "shell.execute_reply": "2024-04-29T14:27:16.463294Z",
     "shell.execute_reply.started": "2024-04-29T14:27:16.438165Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Set up some plotting defaults\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",
    "\n",
    "plt.rcParams.update(params)\n",
    "\n",
    "# Set afw display backend to matplotlib\n",
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "\n",
    "# Set the maximum number of rows to display from pandas\n",
    "pandas.set_option('display.max_rows', 20)\n",
    "\n",
    "plt.style.use('tableau-colorblind10')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad0d6bef-485d-43bd-9ea2-cc79d5579cca",
   "metadata": {},
   "source": [
    "#### 1.2.2 Create some environment variables\n",
    "\n",
    "This next cell creates a new temporary folder in the home directory called \"temp\". The cutout tool must first save cutouts as fits files locally on the Rubin Science Platform (RSP) so this notebook will save those files in temp. At the end of the notebook, the last cell will clear the contents and remove the temp folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "155fd661-08f8-4923-8376-5d6fe5334f4e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:16.464965Z",
     "iopub.status.busy": "2024-04-29T14:27:16.464776Z",
     "iopub.status.idle": "2024-04-29T14:27:16.773065Z",
     "shell.execute_reply": "2024-04-29T14:27:16.772003Z",
     "shell.execute_reply.started": "2024-04-29T14:27:16.464949Z"
    }
   },
   "outputs": [],
   "source": [
    "homeDir = os.getenv(\"HOME\")\n",
    "tempdir = homeDir+'/temp'\n",
    "if not os.path.exists(tempdir):\n",
    "    os.makedirs(tempdir)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff95dea9-1223-4573-adce-bf85285dd2a0",
   "metadata": {},
   "source": [
    "#### 1.2.3 Define functions for manipulating image cutouts and plotting\n",
    "\n",
    "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 1. 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 (cutout service? cutout-sync?)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2f464394-534d-4164-93d4-eaa27f46bae5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:16.777627Z",
     "iopub.status.busy": "2024-04-29T14:27:16.777078Z",
     "iopub.status.idle": "2024-04-29T14:27:17.020521Z",
     "shell.execute_reply": "2024-04-29T14:27:17.019290Z",
     "shell.execute_reply.started": "2024-04-29T14:27:16.777583Z"
    },
    "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()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "16b8a3f9-c6fb-404d-a07c-c852e6afc645",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:17.022977Z",
     "iopub.status.busy": "2024-04-29T14:27:17.022326Z",
     "iopub.status.idle": "2024-04-29T14:27:17.689100Z",
     "shell.execute_reply": "2024-04-29T14:27:17.687938Z",
     "shell.execute_reply.started": "2024-04-29T14:27:17.022930Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def make_image_cutout(tap_service, ra, dec, cutout_size=0.01,\n",
    "                      imtype=None, dataId=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",
    "    cutout_size : 'float'\n",
    "        edge length in degrees of the cutout\n",
    "    imtype : 'string'\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",
    "    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",
    "    filename : 'string'\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/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, s_ra, s_dec \" + \\\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, 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(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",
    "    if filename:\n",
    "        sodaCutout = os.path.join(os.getenv('HOME'), 'temp/'+filename)\n",
    "    else:\n",
    "        sodaCutout = os.path.join(os.getenv('HOME'), 'temp/soda-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": 6,
   "id": "f27addf2-cb60-4456-a6f4-8b8d0bc1b74e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:17.691438Z",
     "iopub.status.busy": "2024-04-29T14:27:17.690812Z",
     "iopub.status.idle": "2024-04-29T14:27:17.959152Z",
     "shell.execute_reply": "2024-04-29T14:27:17.957794Z",
     "shell.execute_reply.started": "2024-04-29T14:27:17.691390Z"
    }
   },
   "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 object from a deepCoadd. \n",
    "\n",
    "\n",
    "### 2.1 Initiate Butler and TAP, and define sky coordinates for the image cutout\n",
    "\n",
    "Both the Butler and the TAP service are needed. The Butler is needed to identify which deepCoadd image hosts the object of interest. The TAP service is then used to query the ivoa.Obscore table for the datalink (a web URL identifying where the data is hosted)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f632c1bc-44dd-4dff-9f29-3e63bf89428d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:17.961405Z",
     "iopub.status.busy": "2024-04-29T14:27:17.960802Z",
     "iopub.status.idle": "2024-04-29T14:27:18.697038Z",
     "shell.execute_reply": "2024-04-29T14:27:18.695980Z",
     "shell.execute_reply.started": "2024-04-29T14:27:17.961360Z"
    }
   },
   "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, we will create a SpherePoint class to define the location on the sky, that is then used to identify which Patch and Tract contains that location among the deepCoadd images.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cdc11b12-7c0b-4ab8-88b0-06787ffd9fe8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-29T14:27:18.699156Z",
     "iopub.status.busy": "2024-04-29T14:27:18.698616Z",
     "iopub.status.idle": "2024-04-29T14:27:19.726838Z",
     "shell.execute_reply": "2024-04-29T14:27:19.724981Z",
     "shell.execute_reply.started": "2024-04-29T14:27:18.699114Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "ename": "DALQueryError",
     "evalue": "Table [ ivoa.ObsCore ] is not found in TapSchema. Possible reasons: table does not exist or permission is denied.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mDALQueryError\u001b[0m                             Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[8], line 17\u001b[0m\n\u001b[1;32m      6\u001b[0m spherePoint \u001b[38;5;241m=\u001b[39m lsst\u001b[38;5;241m.\u001b[39mgeom\u001b[38;5;241m.\u001b[39mSpherePoint(ra\u001b[38;5;241m*\u001b[39mgeom\u001b[38;5;241m.\u001b[39mdegrees, dec\u001b[38;5;241m*\u001b[39mgeom\u001b[38;5;241m.\u001b[39mdegrees)\n\u001b[1;32m      8\u001b[0m query \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSELECT s_ra, s_dec, lsst_patch, lsst_tract, s_region \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \\\n\u001b[1;32m      9\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFROM ivoa.ObsCore \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \\\n\u001b[1;32m     10\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mWHERE dataproduct_type = \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mimage\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \\\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     14\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m, \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mstr\u001b[39m(coord\u001b[38;5;241m.\u001b[39mdec\u001b[38;5;241m.\u001b[39mvalue) \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m), \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \\\n\u001b[1;32m     15\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124ms_region) = 1\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m---> 17\u001b[0m results \u001b[38;5;241m=\u001b[39m \u001b[43mservice\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msearch\u001b[49m\u001b[43m(\u001b[49m\u001b[43mquery\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     19\u001b[0m tract \u001b[38;5;241m=\u001b[39m results[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlsst_tract\u001b[39m\u001b[38;5;124m'\u001b[39m][\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m     20\u001b[0m patch \u001b[38;5;241m=\u001b[39m results[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlsst_patch\u001b[39m\u001b[38;5;124m'\u001b[39m][\u001b[38;5;241m0\u001b[39m]\n",
      "File \u001b[0;32m/opt/lsst/software/stack/conda/miniconda3-py38_4.9.2/envs/lsst-scipipe-8.0.0/lib/python3.11/site-packages/pyvo/dal/tap.py:278\u001b[0m, in \u001b[0;36mTAPService.run_sync\u001b[0;34m(self, query, language, maxrec, uploads, **keywords)\u001b[0m\n\u001b[1;32m    249\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mrun_sync\u001b[39m(\n\u001b[1;32m    250\u001b[0m         \u001b[38;5;28mself\u001b[39m, query, language\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mADQL\u001b[39m\u001b[38;5;124m\"\u001b[39m, maxrec\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, uploads\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m    251\u001b[0m         \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkeywords):\n\u001b[1;32m    252\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    253\u001b[0m \u001b[38;5;124;03m    runs sync query and returns its result\u001b[39;00m\n\u001b[1;32m    254\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    274\u001b[0m \u001b[38;5;124;03m    TAPResults\u001b[39;00m\n\u001b[1;32m    275\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[1;32m    276\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcreate_query\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    277\u001b[0m \u001b[43m        \u001b[49m\u001b[43mquery\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlanguage\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlanguage\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmaxrec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmaxrec\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43muploads\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43muploads\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m--> 278\u001b[0m \u001b[43m        \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkeywords\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mexecute\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/lsst/software/stack/conda/miniconda3-py38_4.9.2/envs/lsst-scipipe-8.0.0/lib/python3.11/site-packages/pyvo/dal/tap.py:1115\u001b[0m, in \u001b[0;36mTAPQuery.execute\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m   1101\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mexecute\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m   1102\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1103\u001b[0m \u001b[38;5;124;03m    submit the query and return the results as a TAPResults instance\u001b[39;00m\n\u001b[1;32m   1104\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1113\u001b[0m \u001b[38;5;124;03m       for errors parsing the VOTable response\u001b[39;00m\n\u001b[1;32m   1114\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1115\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mTAPResults\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mexecute_votable\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43murl\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mqueryurl\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msession\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_session\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/lsst/software/stack/conda/miniconda3-py38_4.9.2/envs/lsst-scipipe-8.0.0/lib/python3.11/site-packages/pyvo/dal/adhoc.py:111\u001b[0m, in \u001b[0;36mAdhocServiceResultsMixin.__init__\u001b[0;34m(self, votable, url, session)\u001b[0m\n\u001b[1;32m    110\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;28mself\u001b[39m, votable, url\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, session\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 111\u001b[0m     \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mvotable\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43murl\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43murl\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msession\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msession\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    113\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_adhocservices \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(\n\u001b[1;32m    114\u001b[0m         resource \u001b[38;5;28;01mfor\u001b[39;00m resource \u001b[38;5;129;01min\u001b[39;00m votable\u001b[38;5;241m.\u001b[39mresources\n\u001b[1;32m    115\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m resource\u001b[38;5;241m.\u001b[39mtype \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmeta\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m resource\u001b[38;5;241m.\u001b[39mutype \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124madhoc:service\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    116\u001b[0m     )\n",
      "File \u001b[0;32m/opt/lsst/software/stack/conda/miniconda3-py38_4.9.2/envs/lsst-scipipe-8.0.0/lib/python3.11/site-packages/pyvo/dal/query.py:322\u001b[0m, in \u001b[0;36mDALResults.__init__\u001b[0;34m(self, votable, url, session)\u001b[0m\n\u001b[1;32m    320\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_status \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_findstatus(votable)\n\u001b[1;32m    321\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_status[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mlower() \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m (\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mok\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124moverflow\u001b[39m\u001b[38;5;124m\"\u001b[39m):\n\u001b[0;32m--> 322\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m DALQueryError(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_status[\u001b[38;5;241m1\u001b[39m], \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_status[\u001b[38;5;241m0\u001b[39m], url)\n\u001b[1;32m    324\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_status[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mlower() \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124moverflow\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m    325\u001b[0m     warn(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPartial result set. Potential causes MAXREC, async storage space, etc.\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    326\u001b[0m          category\u001b[38;5;241m=\u001b[39mDALOverflowWarning)\n",
      "\u001b[0;31mDALQueryError\u001b[0m: Table [ ivoa.ObsCore ] is not found in TapSchema. Possible reasons: table does not exist or permission is denied."
     ]
    }
   ],
   "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",
    "\n",
    "query = \"SELECT s_ra, s_dec, lsst_patch, lsst_tract, s_region \" + \\\n",
    "        \"FROM ivoa.ObsCore \" + \\\n",
    "        \"WHERE dataproduct_type = 'image' \" + \\\n",
    "        \"AND obs_collection = 'LSST.DP02' \" + \\\n",
    "        \"AND dataproduct_subtype = 'lsst.deepCoadd_calexp' \" + \\\n",
    "        \"AND CONTAINS(POINT('ICRS', \" + str(coord.ra.value) + \\\n",
    "        \", \" + str(coord.dec.value) + \"), \" + \\\n",
    "        \"s_region) = 1\"\n",
    "\n",
    "results = service.search(query)\n",
    "\n",
    "tract = results['lsst_tract'][0]\n",
    "patch = results['lsst_patch'][0]\n",
    "\n",
    "# check datalab tutorials for more info about SODA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "248c8160-072b-4497-8820-0e979162fd8d",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f075025-e609-4c1c-a965-eda296fcd549",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.727466Z",
     "iopub.status.idle": "2024-04-29T14:27:19.727722Z",
     "shell.execute_reply": "2024-04-29T14:27:19.727603Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.727591Z"
    }
   },
   "outputs": [],
   "source": [
    "dataId = {'band': 'i', 'tract': tract,\n",
    "          'patch': patch}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "551c859a-730a-4f4e-811b-6c2a13c4c30d",
   "metadata": {},
   "source": [
    "The next cell shows the query to the TAP service 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 images from observations is to use the TAP service to query the ObsCore schema. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2991ef70-d628-4c2f-b2d4-5f403a798556",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.728949Z",
     "iopub.status.idle": "2024-04-29T14:27:19.729187Z",
     "shell.execute_reply": "2024-04-29T14:27:19.729078Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.729068Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "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(tract) + \" \" + \\\n",
    "    \"AND lsst_patch = \" + str(patch) + \" \" + \\\n",
    "    \"AND lsst_band = 'i' \"\n",
    "\n",
    "results = service.search(query)\n",
    "results.to_table().show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "996dd37e-7c19-4796-8b19-5362d347bc29",
   "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",
    "### 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 (see DP0.2 Tutorial 22a, which has further technical explanation). 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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.729886Z",
     "iopub.status.idle": "2024-04-29T14:27:19.730103Z",
     "shell.execute_reply": "2024-04-29T14:27:19.730005Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.729996Z"
    },
    "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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.731396Z",
     "iopub.status.idle": "2024-04-29T14:27:19.731633Z",
     "shell.execute_reply": "2024-04-29T14:27:19.731528Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.731519Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = SodaQuery.from_resource(dl_results,\n",
    "                             dl_results.get_adhocservice_by_id(\"cutout-sync\"),\n",
    "                             session=auth_session)\n"
   ]
  },
  {
   "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 2 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 rectanglar 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. \n",
    "\n",
    "The output of the cutout service is (currently) limited to a fits file stored in the returned Google cloud storage service (GCS) bucket. Thus, while the image cutout procedure itself is performed remotely, the cutout must be written to disk first before it is accessed by the notebook (which is done using the write() function). (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 will be stored in the $HOME/temp folder. Finally, call the plotImage function (defined above), which will quickly display the cutout that was written to disk in in the temp folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23ca3f11-f17e-4e5c-8469-a0d25554234d",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.732281Z",
     "iopub.status.idle": "2024-04-29T14:27:19.732557Z",
     "shell.execute_reply": "2024-04-29T14:27:19.732443Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.732434Z"
    }
   },
   "outputs": [],
   "source": [
    "sphereRadius = 0.05 * u.deg\n",
    "\n",
    "sq.circle = (spherePoint.getRa().asDegrees() * u.deg,\n",
    "             spherePoint.getDec().asDegrees() * u.deg,\n",
    "             sphereRadius)\n",
    "\n",
    "# first write the file to disk:\n",
    "sodaCutout = os.path.join(os.getenv('HOME'), 'temp/soda-cutout-circle.fits')\n",
    "with open(sodaCutout, 'bw') as f:\n",
    "    f.write(sq.execute_stream().read())\n",
    "\n",
    "# Display the cutout from the file\n",
    "plotImage(ExposureF(sodaCutout))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99e35e80-15ce-40c6-aa8b-9079e3bc0200",
   "metadata": {},
   "source": [
    "### 2.2.2 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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.733344Z",
     "iopub.status.idle": "2024-04-29T14:27:19.733581Z",
     "shell.execute_reply": "2024-04-29T14:27:19.733472Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.733462Z"
    },
    "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",
    "#               spherePoint.getRa().asDegrees()* u.deg + sphereRadius2,\n",
    "#               spherePoint.getDec().asDegrees()*u.deg - sphereRadius2)\n",
    "\n",
    "sodaCutout2 = os.path.join(os.getenv('HOME'), 'temp/soda-cutout-polygon.fits')\n",
    "with open(sodaCutout2, 'bw') as f:\n",
    "    f.write(sq2.execute_stream().read())\n",
    "\n",
    "plotImage(ExposureF(sodaCutout2))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6970143e-90b9-406c-80ce-e6c6ba923c1e",
   "metadata": {},
   "source": [
    "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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.734679Z",
     "iopub.status.idle": "2024-04-29T14:27:19.734921Z",
     "shell.execute_reply": "2024-04-29T14:27:19.734805Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.734795Z"
    }
   },
   "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'), 'temp/soda-cutout_polygon.fits')\n",
    "with open(sodaCutout2, 'bw') as f:\n",
    "    f.write(sq2.execute_stream().read())\n",
    "\n",
    "plotImage(ExposureF(sodaCutout2))\n"
   ]
  },
  {
   "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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.735745Z",
     "iopub.status.idle": "2024-04-29T14:27:19.735987Z",
     "shell.execute_reply": "2024-04-29T14:27:19.735887Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.735877Z"
    },
    "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))\n"
   ]
  },
  {
   "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. 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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.736781Z",
     "iopub.status.idle": "2024-04-29T14:27:19.737020Z",
     "shell.execute_reply": "2024-04-29T14:27:19.736916Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.736906Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "center_coords = \"62, -37\"\n",
    "ralbg = 62.\n",
    "declbg = -37.\n",
    "radius = \"4.0\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c3e8bae-e010-4b2e-b3fc-e688ebf42c5b",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.737962Z",
     "iopub.status.idle": "2024-04-29T14:27:19.738190Z",
     "shell.execute_reply": "2024-04-29T14:27:19.738085Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.738076Z"
    }
   },
   "outputs": [],
   "source": [
    "max_rec = 500000\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",
    "\n",
    "results = service.search(query)\n",
    "\n",
    "results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77189c9e-85b8-441f-9d93-ba1277bb598d",
   "metadata": {},
   "source": [
    "The next cell will plot 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).\n",
    "\n",
    "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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe52da8d-08e7-43a8-83ea-ccfe32e1ec10",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.739144Z",
     "iopub.status.idle": "2024-04-29T14:27:19.739382Z",
     "shell.execute_reply": "2024-04-29T14:27:19.739275Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.739265Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#  here, set u-band magnitude to umag of 1-sigma error before selecting colors\n",
    "umag = np.zeros(len(results['umag']))\n",
    "ubandS2N = results['u_ap09Flux'] / results['u_ap09FluxErr']\n",
    "\n",
    "for i in range(len(umag)):\n",
    "    whnondet = np.where(ubandS2N[i] < 1)[0]\n",
    "    if len(whnondet) > 0:\n",
    "        umag[i] = results['umagErr'][i]\n",
    "    else:\n",
    "        umag[i] = results['umag'][i]\n",
    "\n",
    "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",
    "# plot the color-color selection box:\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()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f93f4a8a-d624-4d68-8591-0ebe4d0c5a8f",
   "metadata": {},
   "source": [
    "Some fraction of the orange LBG candidates might be interlopers. The cell below will plot image cutouts of each LBG candidate for visual inspection. Since this is a simulated dataset, included in this cell is `query3` to the TAP service, which pulls the intrinsic redshift for each object from the truth table to compare with the expected redshift $z\\sim3$ of selected objects (see DP0.2 Tutorial Notebook 08 for more information on truth tables). \n",
    "\n",
    "Due to the number of cutouts, this cell can take some time to run. Adjust the range of the for loop as needed to see the cutouts for the entire sample, rather than a subset of 5 (for example, to display the full set of selected LBGs, replace the first line with `for i in range(len(whlbg)):`).\n",
    "\n",
    "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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.740349Z",
     "iopub.status.idle": "2024-04-29T14:27:19.740577Z",
     "shell.execute_reply": "2024-04-29T14:27:19.740474Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.740464Z"
    }
   },
   "outputs": [],
   "source": [
    "filt = [\"u\", \"g\", \"r\", \"i\", \"z\"]\n",
    "\n",
    "for i in range(5):\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",
    "        id = str(results['objectId'][whlbg][i])\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=\"+id+\" \" + \\\n",
    "                 \"AND ts.truth_type=1 \" + \\\n",
    "                 \"AND obj.detect_isPrimary=1 \"\n",
    "\n",
    "        results3 = service.search(query3)\n",
    "        results3.to_table()\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'][0]),\n",
    "                color='white')\n"
   ]
  },
  {
   "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": "cc15bed2-16e3-421d-8f52-798f01090973",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.741539Z",
     "iopub.status.idle": "2024-04-29T14:27:19.741760Z",
     "shell.execute_reply": "2024-04-29T14:27:19.741665Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.741656Z"
    }
   },
   "outputs": [],
   "source": [
    "diaObjectId = 1250953961339360185\n",
    "\n",
    "results = service.search(\"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",
    "\n",
    "diaSrc = results.to_table()\n",
    "diaSrc"
   ]
  },
  {
   "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": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.742488Z",
     "iopub.status.idle": "2024-04-29T14:27:19.742720Z",
     "shell.execute_reply": "2024-04-29T14:27:19.742612Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.742603Z"
    },
    "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_fitsname = 'cutout_'+str(i)+'.fits'\n",
    "    image_file = os.path.join(os.getenv('HOME'), 'temp/'+ref_fitsname)\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\n"
   ]
  },
  {
   "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 eachother, 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 eachother. 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). (Note: ExposureF creates a known warning that is benign in this circumstance; it originates because the recommended LSST science pipeline version used here, 2024_4, is behind the version used when creating cutouts which is 2024_5. The warning can safely be ignored).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88303981-d84d-45f2-b1ad-6a6a761b3141",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.743980Z",
     "iopub.status.idle": "2024-04-29T14:27:19.744212Z",
     "shell.execute_reply": "2024-04-29T14:27:19.744105Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.744096Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "i = 0\n",
    "cutout_size_pix = 131\n",
    "\n",
    "# 0.2-arcsec pixels, convert to degrees\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_fitsname = 'cutout_ref.fits'\n",
    "        ref_cutout = os.path.join(os.getenv('HOME'), 'temp/'+ref_fitsname)\n",
    "\n",
    "        make_image_cutout(service, ra, dec, cutout_size=cutout_size,\n",
    "                          imtype='calexp', dataId=dataId_calexp,\n",
    "                          filename=ref_fitsname)\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",
    "                               'temp/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'), 'temp/'+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",
    "                               'temp/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\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b921539-cd37-4cd1-95d4-fb1a5ffca9e2",
   "metadata": {},
   "source": [
    "### 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. Once that is completed, the second cell below will do the file cleanup, and delete all of the fits, png, and gifs created during this tutorial. To remove the temp folder, close the notebook, open the terminal and cd to the home directory, and type `rmdir temp` from the command line.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78a1cb9e-c0cb-46e3-af34-11360b8a6b1a",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.745040Z",
     "iopub.status.idle": "2024-04-29T14:27:19.745264Z",
     "shell.execute_reply": "2024-04-29T14:27:19.745160Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.745151Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "images = os.path.join(os.getenv('HOME'), 'temp/')\n",
    "make_gif(images)\n",
    "\n",
    "display(dimg(data=open(images+'animation.gif', 'rb').read()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00d1d56b-3dbd-4bb9-8799-771f474c2725",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-04-29T14:27:19.745950Z",
     "iopub.status.idle": "2024-04-29T14:27:19.746182Z",
     "shell.execute_reply": "2024-04-29T14:27:19.746076Z",
     "shell.execute_reply.started": "2024-04-29T14:27:19.746066Z"
    }
   },
   "outputs": [],
   "source": [
    "for filepath in glob.glob(os.path.join(os.getenv('HOME'), 'temp/')+\"*\"):\n",
    "    os.remove(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0138cbf3-11d9-465e-9b62-7de23ed2e165",
   "metadata": {},
   "source": [
    "## 5. 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.7"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}