cutout_example.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>Image Cutouts</b> <br>\n",
    "Contact author(s): Christina Williams <br>\n",
    "Last verified to run: 2023-09-12 <br>\n",
    "LSST Science Piplines version: Weekly 2023_21 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** Demonstration of how to use the cutout tool for calexp and deepCoadds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This is a brief demonstration of how to make a single-filter image cutout using the cutout tool.\n",
    "\n",
    "The calls to the cutout tool are embedded in the `makeImageCutout` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1 Package Imports\n",
    "\n",
    "Import general python and astronomy packages, LSST Science Pipelines packages, and image cutout `PyVO` packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy import units as u\n",
    "\n",
    "from lsst.daf.butler import Butler\n",
    "import lsst.geom as geom\n",
    "from lsst.afw.image import Image\n",
    "from lsst.afw.image import Exposure, ExposureF\n",
    "import lsst.afw.display as afwDisplay\n",
    "from lsst.rsp import get_tap_service\n",
    "\n",
    "from pyvo.dal.adhoc import DatalinkResults, SodaQuery"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9b2f7b1-691a-47e2-8efd-292c2b059d38",
   "metadata": {},
   "source": [
    "The default option is that image cutouts are saved as `~/DATA/soda-cutout.fits`.\n",
    "\n",
    "Image cutouts can be renamed and stored in the directory `./cutouts`. If that directory does not exist, create it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f676e80-8191-4147-9dda-5fb03d31cdd8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "if os.path.exists('cutouts/'):\n",
    "    print('Exists: cutouts/')\n",
    "else:\n",
    "    os.system('mkdir cutouts')\n",
    "    print('Created: cutouts/')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80bedd3d-8d87-4eb8-a1df-50c0c478a6a6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-12T20:20:13.890045Z",
     "iopub.status.busy": "2023-09-12T20:20:13.889302Z",
     "iopub.status.idle": "2023-09-12T20:20:13.894641Z",
     "shell.execute_reply": "2023-09-12T20:20:13.893900Z",
     "shell.execute_reply.started": "2023-09-12T20:20:13.890018Z"
    },
    "tags": []
   },
   "source": [
    "Set the `afwDisplay` backend to be `matplotlib`, and the maximum number of rows to display for a `pandas` table to 20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc025155-2064-445d-a7da-f17f80982f04",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "pandas.set_option('display.max_rows', 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d48528ff-f907-4b0b-b8e3-8736b967df07",
   "metadata": {},
   "source": [
    "Function `plotImage` is simply for convenience of plotting the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8503823e-780f-4dbe-9209-d582b3acf538",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotImage(exposure: ExposureF, fig_width=5, fig_height=5):\n",
    "    \"\"\"Plot an image using matplotlib\n",
    "   \n",
    "    Parameters\n",
    "    ----------\n",
    "    image : `Exposure`\n",
    "        the image to plot\n",
    "    fig_width : int\n",
    "        width of display of cutout, does not affect cutout size\n",
    "    fig_height : int\n",
    "        height of display of cutout, does not affect cutout size\n",
    "    \"\"\"\n",
    "    \n",
    "    fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "    display = afwDisplay.Display(frame=fig)\n",
    "    display.scale('asinh', 'zscale')\n",
    "    display.mtv(exposure.image)\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "115935db-38d1-4a92-9551-0ef116188671",
   "metadata": {},
   "source": [
    "Function `makeImageCutout` provides the cutout service for either `calexp` or `deepCoadds` images.\n",
    "\n",
    "Note that `SodaQuery.from_resource` creates a instance from a number of records and a Datalink Resource."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66431945-0740-475d-9256-492f6a7b633c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def makeImageCutout(tap_service, ra, dec, cutout_size=0.01, default_filter = 'i', \n",
    "                    imtype='deepCoadd', dataId=None, filename=None): \n",
    "\n",
    "    \"\"\"Make a cutout using the cutout tool\n",
    "   \n",
    "    Parameters\n",
    "    ----------\n",
    "    tap_service : \n",
    "       the tap service instantiated with get_tap_service(\"tap\")\n",
    "    \n",
    "    ra : `float`\n",
    "       right asencion of the cutout center, in decimal degrees\n",
    "    \n",
    "    dec : `float`\n",
    "       declination of the cutoutcenter, in decimal degrees\n",
    "    \n",
    "    cutout_size : `float`\n",
    "       the size of the side of the cutout in degrees\n",
    "    \n",
    "    default_filter : `str`\n",
    "       the default filter to use in a deepCoadd dataId if it is not supplied\n",
    "    \n",
    "    imtype : `str`\n",
    "       calexp or deepCoadd (default)\n",
    "    \n",
    "    dataId : dict\n",
    "       dataId to pass to the bulter\n",
    "       a deepCoadd dataId should have format {'band':'band', 'tract':tract, 'patch':patch}\n",
    "       a calexp dataId should have format {'visit':visitId, 'detector':detector}\n",
    "    \n",
    "    filename : `str`\n",
    "       to save the cutout as a FITS file in \"cutouts/filename.fits\"\n",
    "       otherwise it is saved in \"~/DATA/soda-cutout.fits\"\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    sodaCutout : `str`\n",
    "       the filename of the saved image cutout\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    spherePoint = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)\n",
    "    \n",
    "    if imtype == 'calexp':\n",
    "        query = \"SELECT access_format, access_url, dataproduct_subtype, lsst_visit, lsst_detector, \" + \\\n",
    "            \"lsst_band, s_ra, s_dec 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",
    "    else:\n",
    "        tract = dataId[\"tract\"]\n",
    "        patch = dataId[\"patch\"]\n",
    "        if 'band' in dataId:\n",
    "            band = dataId[\"band\"]\n",
    "        else:\n",
    "            band = default_filter\n",
    "        query = \"SELECT access_format, access_url, dataproduct_subtype, lsst_patch, lsst_tract, \" + \\\n",
    "            \"lsst_band, s_ra, s_dec 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",
    "    dataLinkUrl = results[0].getdataurl()\n",
    "    f\"Datalink link service url: {dataLinkUrl}\"\n",
    "    auth_session = service._session\n",
    "    dl_results = DatalinkResults.from_result_url(dataLinkUrl,session=auth_session)\n",
    "    f\"{dl_results.status}\"\n",
    "\n",
    "    sq = SodaQuery.from_resource(dl_results, dl_results.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 = 'cutouts/'+filename       \n",
    "    else:\n",
    "        sodaCutout = os.path.join(os.getenv('HOME'), 'DATA/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": "ec51ac0b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Make a `deepCoadd` image cutout\n",
    "\n",
    "Instantiate access to the data butler and the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b8cd59-1ba3-4eaa-846f-6478ed0c3cf5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "butler = Butler('dp02', collections='2.2i/runs/DP0.2')\n",
    "registry = butler.registry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a02ecd1-f6d0-45e5-9f89-e1dc56f16eff",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af82fa77-447e-4dad-a5ee-e0f9c6a677c9",
   "metadata": {
    "tags": []
   },
   "source": [
    "Define the `dataId` for a given sky coordinate (in this case, the location of a known galaxy cluster) and filter (band).\n",
    "\n",
    "Use the butler `skyMap` to obtain the `tract` and `patch` info, which are required for the `dataId`, for these coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc11b12-7c0b-4ab8-88b0-06787ffd9fe8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ra = 55.7467 \n",
    "dec = -32.2862\n",
    "spherePoint = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)\n",
    "skymap = butler.get('skyMap')\n",
    "tract = skymap.findTract(spherePoint)\n",
    "patch = tract.findPatch(spherePoint)\n",
    "dataId = {'band':'i', 'tract':tract.tract_id, 'patch': patch.getSequentialIndex()}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c53ba0bc-91c2-4576-8014-4e008929d50c",
   "metadata": {},
   "source": [
    "Create and display the `deepCoadd` image cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd1a68b4-5524-4016-a743-88b34e1d661a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_cutout = makeImageCutout(service, ra, dec, dataId=dataId, cutout_size=0.1)\n",
    "plotImage(ExposureF(my_cutout))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eac08b14-be96-44eb-80a2-7803b9914412",
   "metadata": {},
   "source": [
    "Find the above cutout saved as a FITS file in `~/DATA/soda-cutout.fits`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5111d5-82a2-43ec-b4a8-daafdda10e72",
   "metadata": {},
   "source": [
    "## 3. Make a series of `calexp` image cutouts\n",
    "\n",
    "In this example, for a given `DiaObject`, create small image cutouts (stamps) from the first 10 _i_-band `calexp` images ever obtained by LSST.\n",
    "\n",
    "This example does not limit cutout creation only to the visits in which the `DiaObject` was detected in the difference-image, so there might not be a source at the cutout's center.\n",
    "\n",
    "Use the TAP service to obtain the `visitId` and `detector` identifiers for images containing this `DiaObject`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bdf25f3-f97d-44fa-91e3-8f66f18d6953",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diaObjectId = 1253478440036730088\n",
    "\n",
    "ccdquery = \"SELECT TOP 10 dia.coord_ra, dia.coord_dec, \" + \\\n",
    "           \"dia.diaObjectId,  dia.ccdVisitId, dia.band, \" + \\\n",
    "           \"cv.visitId, cv.physical_filter, cv.detector, cv.obsStartMJD, \" + \\\n",
    "           \"cv.expMidptMJD \" + \\\n",
    "           \"FROM dp02_dc2_catalogs.ForcedSourceOnDiaObject as dia \" + \\\n",
    "           \"JOIN dp02_dc2_catalogs.CcdVisit as cv ON cv.ccdVisitId = dia.ccdVisitId \" + \\\n",
    "           \"WHERE dia.diaObjectId = \"+str(diaObjectId)+\" AND cv.band = 'i'\"\n",
    "    \n",
    "diaccdsearch = service.search(ccdquery)\n",
    "forcedSrc = diaccdsearch.to_table()\n",
    "del diaccdsearch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edc4ab68-01ae-408e-a5f4-01299d0bb3b4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-12T20:36:30.026103Z",
     "iopub.status.busy": "2023-09-12T20:36:30.025276Z",
     "iopub.status.idle": "2023-09-12T20:36:30.030582Z",
     "shell.execute_reply": "2023-09-12T20:36:30.029843Z",
     "shell.execute_reply.started": "2023-09-12T20:36:30.026072Z"
    },
    "tags": []
   },
   "source": [
    "Sort the visits by MJD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05226e70-d21e-4a2b-972b-2ffc09daa1d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wh = np.argsort(forcedSrc['obsStartMJD'])\n",
    "forcedSrc[wh]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29bf7996-7ffb-4a5e-9367-d899180880f3",
   "metadata": {},
   "source": [
    "Create 10 small cutouts, display them, and save them with filenames `cutouts/cutout_N.fits` where N will be 0 to 9.\n",
    "\n",
    "Note: the cutouts are not uniformly North up East left."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a5bdf81-8438-4cb8-806b-d4e0778da92c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in range(10):\n",
    "    \n",
    "    ra = forcedSrc['coord_ra'][wh][i]\n",
    "    dec = forcedSrc['coord_dec'][wh][i]\n",
    "    spherePoint = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)\n",
    "    dataId_calexp = {'visit':forcedSrc['visitId'][wh][i], 'detector':forcedSrc['detector'][wh][i]}\n",
    "\n",
    "    my_cutout = makeImageCutout(service, ra, dec, cutout_size=0.002, imtype='calexp', dataId=dataId_calexp,\n",
    "                                filename='cutout_'+str(i)+'.fits')\n",
    "    \n",
    "    plotImage(ExposureF(my_cutout), fig_width=2, fig_height=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "841f517b-2480-418d-9829-ad2b0c79d148",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}