DP02_03c_Big_deepCoadd_Cutout.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "973d2e4b-1fb6-4af3-be5e-f04fd1b3390a",
   "metadata": {},
   "source": [
    "# Create a large cutout of coadded images\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 authors: Christina Williams, Melissa Graham <br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_50 <br>\n",
    "Container Size: Large <br>\n",
    "Targeted learning level: Intermediate <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc18d3c3-0a18-4e6f-9f5e-ec4b0bb7dc64",
   "metadata": {},
   "source": [
    "**Description:** Create a large custom `deepCoadd` cutout."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de444327-7b83-48dd-af78-1b767f30e92c",
   "metadata": {},
   "source": [
    "**Skills:** Identify tracts and patches and combine them into a large custom `deepCoadd` cutout."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65848064-fe12-4060-8650-3ad786309402",
   "metadata": {},
   "source": [
    "**LSST Data Products:** deepCoadd_calexp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14a042a2-5152-4256-9190-7a81222e1c1f",
   "metadata": {},
   "source": [
    "**Packages:** lsst.ip.diffim.GetTemplateTask, lsst.afw.geom.makeSkyWcs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce2ecd7d-2b7d-4fea-8501-7b9a6f6b5df3",
   "metadata": {},
   "source": [
    "**Credit:** Thanks to Aline Chu for raising in the Rubin Community Forum the issue of how to create patch- and tract-spanning cutouts. Thanks to Gregory Dubois-Felsmann and Lauren MacArthur who helped in the Forum, and John Parejko, Yusra Al-Sayyad, Tim Jenness and others for help with the LSST science pipelines."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5c4f690-88d7-4d69-aea5-0bc255372db8",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>.\n",
    "Questions are welcome as new topics in the \n",
    "<a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> \n",
    "of the Rubin Community Forum. \n",
    "Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89c56237-1162-40f3-bac4-0171f6f61944",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "At the moment, the image cutout service demonstrated in DP0.2 tutorial notebook\n",
    "13a titled <a href=\"https://github.com/rubin-dp0/tutorial-notebooks/blob/main/DP02_13a_Image_Cutout_SciDemo.ipynb\">DP02_13a_Image_Cutout_SciDemo.ipynb</a> must be passed a tract and patch, and can only create a cutout from a single\n",
    "patch.\n",
    "\n",
    "This will always be a contraint for single processed visit images (PVIs,\n",
    "also known as `calexp`s).\n",
    "In the future, users will be able to request `deepCoadd` cutouts from the \n",
    "Rubin image cutout service which span patch and tract boundaries.\n",
    "\n",
    "In the meantime, this tutorial demonstrates how to use the `GetTemplateTask`\n",
    "to create a single cutout image with contributions from multiple adjacent\n",
    "patches and tracts. \n",
    "\n",
    "In the LSST Science Pipelines, `GetTemplateTask` is used to create a template\n",
    "image from `deepCoadd` images for a given PVI, in order to perform Difference Image Analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44158545-a490-4c7a-a92e-ed951f943db1",
   "metadata": {},
   "source": [
    "These links provide documentation for\n",
    "<a href=\"https://pipelines.lsst.io/modules/lsst.ip.diffim/tasks/lsst.ip.diffim.GetTemplateTask.html\">GetTemplateTask</a>\n",
    "and its\n",
    "<a href=\"https://pipelines.lsst.io/py-api/lsst.ip.diffim.GetTemplateTask.html#lsst.ip.diffim.GetTemplateTask.run\">run method</a>. To summarize,\n",
    "the function of this task is to build a \"template\" combined image from existing coadds. In the simplest case where the coadds that contribute to the assembled template image are all patches of the same tract, the WCS and bounding box of the pixel map is the same across the patches. The documentation for `GetTemplateTask` for this simple case of joining patches from a common tract says:\n",
    ">_\"...overlapping box regions of the input template patches are pixel by pixel copied into the assembled template image. There is no warping or pixel resampling. Pixels with no overlap of any available input patches are set to nan value and NO_DATA flagged.\"_\n",
    "\n",
    "However, for the more general case where input coadds may come from different tracts, the input images have different WCS. Thus, `GetTemplateTask` will warp the coadds from different tracts onto a common projection (i.e. a newly defined WCS). This tutorial will also demonstrate how to generate the new common WCS using `makeSkyWcs`.  \n",
    "\n",
    "The `GetTemplateTask` documentation further notes:\n",
    ">_\"Where the tracts overlap, the resulting template image is averaged. The PSF on the template is created by combining the CoaddPsf on each template image into a meta-CoaddPsf.\"_\n",
    "\n",
    "This notebook will demonstrate how to access the final data products including the `CoaddPsf` in Section 5.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb6377f3-8a35-451e-9fe3-a8989a944f53",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-04T21:56:55.196119Z",
     "iopub.status.busy": "2024-05-04T21:56:55.195330Z",
     "iopub.status.idle": "2024-05-04T21:56:55.199978Z",
     "shell.execute_reply": "2024-05-04T21:56:55.199425Z",
     "shell.execute_reply.started": "2024-05-04T21:56:55.196089Z"
    }
   },
   "source": [
    "**Related tutorials:**\n",
    "This tutorial assumes some familiarity with the following concepts.\n",
    "DP0.2 tutorial notebooks that demonstrate these concepts more fully are listed in parentheses.\n",
    "It is not necessary to run through all of these other tutorials before doing this one,\n",
    "but they might be useful as reference.\n",
    "\n",
    " * the butler (04a, 04b)\n",
    " * the TAP service (02)\n",
    " * image display (03a, 03b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a3f0133-53e2-4fb4-9643-8d8351b9ef26",
   "metadata": {},
   "source": [
    "### 1.1. Import packages\n",
    "\n",
    "Import packages from the LSST Science Pipelines (top) and \n",
    "general python packages (bottom)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0195d52-71b0-4d34-b7af-9988cc522a26",
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.geom as geom\n",
    "from lsst.rsp import get_tap_service, retrieve_query\n",
    "from lsst.daf.butler import Butler, CollectionType\n",
    "from lsst.ip.diffim import GetTemplateTask\n",
    "from lsst.afw.image import Image\n",
    "from lsst.afw.geom import makeSkyWcs\n",
    "\n",
    "import collections\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from astropy.wcs import WCS\n",
    "\n",
    "from photutils.aperture import SkyCircularAperture\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "456ae0ad-2f86-4251-9e24-fcaabf32af18",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters\n",
    "\n",
    "Set the `afwDisplay` backend to `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "618960a9-c365-438c-bd25-737fdc8f9055",
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e29a746-6a8b-4087-b21d-96949be20661",
   "metadata": {},
   "source": [
    "Instantiate the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "974cae62-db64-469b-91c2-1fcb2b032949",
   "metadata": {},
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47c63680-cf3d-485c-a2d1-fadbf1e2c296",
   "metadata": {},
   "source": [
    "Instantiate the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da77a370-98df-44fa-903a-a26b69e83235",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = 'dp02' \n",
    "collection = '2.2i/runs/DP0.2'\n",
    "butler = Butler(config, collections=collection)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f00a85e-c64f-473b-84e7-df0feaa0bdfc",
   "metadata": {},
   "source": [
    "## 2. Prepare to create the cutout\n",
    "\n",
    "### 2.1. Define central coordinates and box size\n",
    "\n",
    "Define the central coordinates and the size of the side of the desired cutout.\n",
    "\n",
    "For this example, the coordinates 50.5265, -39.7589 are a galaxy cluster known to be near the intersection of patches and tracts, but whose extent exceeds the distance to the image edges.\n",
    "\n",
    "The size of the cutout will be constrained by the container size selected\n",
    "when the Notebook Aspect session was started, because all the images\n",
    "that contribute to the cutout must be held in memory simultaneously\n",
    "and passed to `GetTemplateTask` (this will become clear in Section 3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4918ac1-7798-494e-b78d-d8b0e16d5b6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_ra = 50.5265\n",
    "my_dec = -39.7589\n",
    "box_side_deg = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ed96627-df56-4e30-8d22-747a81ba0fb6",
   "metadata": {},
   "source": [
    "Define the extent of the cutout in RA and Dec. Include the `cos(dec)` factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c3239f8-b753-42f4-b43a-d342f14c1fe0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ra1 = my_ra + 0.5 * box_side_deg / np.cos(np.radians(my_dec))\n",
    "ra2 = my_ra - 0.5 * box_side_deg / np.cos(np.radians(my_dec))\n",
    "dec1 = my_dec - 0.5 * box_side_deg\n",
    "dec2 = my_dec + 0.5 * box_side_deg\n",
    "print('ra1 = %8.4f,  ra2 = %8.4f,  delta RA = %6.3f' %\n",
    "      (ra1, ra2, np.abs(ra2-ra1)))\n",
    "print('dec1 = %8.4f, dec2 = %8.4f, delta Dec = %6.3f' %\n",
    "      (dec1, dec2, np.abs(dec2-dec1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b051dd23-a6b2-480c-841b-d06b1f4cce93",
   "metadata": {},
   "source": [
    "Define corners of the cutout.\n",
    "\n",
    "The order is: SE, SW, NW, NE, as drawn in the plot below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0017729-c154-44a4-b005-168072327cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "corners = np.asarray([[ra1, dec1], [ra2, dec1],\n",
    "                      [ra2, dec2], [ra1, dec2]],\n",
    "                     dtype='float')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "437d51e3-8b17-40c1-8488-95ec42d17cac",
   "metadata": {},
   "source": [
    "Create a little plot to visualize the corners.\n",
    "\n",
    "Note that the figure is forced to 3:2 dimensions size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5df7472-f144-4011-a9fc-dda77a35815c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3, 2))\n",
    "plt.plot(my_ra, my_dec, '*', ms=10, color='blue')\n",
    "clabels = ['1. SE', '2. SW', '3. NW', '4. NE']\n",
    "for corner, clabel in zip(corners, clabels):\n",
    "    plt.text(corner[0], corner[1], clabel, ha='center', va='center')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.xlim([ra1+0.05, ra2-0.05])\n",
    "plt.ylim([dec1-0.05, dec2+0.05])\n",
    "plt.xlabel('RA')\n",
    "plt.ylabel('Dec')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23ec8b4b-948a-49ae-830c-3dc37246f641",
   "metadata": {},
   "source": [
    "> Figure 1: A small plot of RA vs. Dec to visualize the center (blue star) and the four corners (numbered and labeled) defined by the user for the large cutout. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63e545bc-1c79-4b38-afe4-e8fc357e5e2b",
   "metadata": {},
   "source": [
    "### 2.2. Identify the patch to use as anchor image\n",
    "\n",
    "Use the butler to find the tract and patch within which\n",
    "the central coordinates of the desired cutouts are most centered.\n",
    "\n",
    "This will be the tract and patch referred to as the \"anchor image\"\n",
    "when creating the big cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2246af6-08dd-48fa-8a1c-cdc3eb17449b",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_point = geom.SpherePoint(my_ra * geom.degrees,\n",
    "                            my_dec * geom.degrees)\n",
    "skymap = butler.get('skyMap')\n",
    "tract = skymap.findTract(my_point)\n",
    "patch = tract.findPatch(my_point)\n",
    "anchor_tract = tract.tract_id\n",
    "anchor_patch = patch.getSequentialIndex()\n",
    "print('anchor tract and patch: ', anchor_tract, anchor_patch)\n",
    "del my_point, skymap, tract, patch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cee06e4-087a-4fcb-b291-e6f16bda8c90",
   "metadata": {},
   "source": [
    "### 2.3. Identify additional tracts and patches\n",
    "\n",
    "In the future, it will be possible to write ADQL statements\n",
    "that use the INTERSECTS functionality to identify all patches\n",
    "(or images) that overlap with a defined polygon.\n",
    "\n",
    "For now, use the TAP service to search the `coaddPatches` table\n",
    "for images with central coordinates that indicate the patch\n",
    "could *potentially* overlap with the defined extent of the \n",
    "desired cutout.\n",
    "\n",
    "Increase the extent of the search by the patch size to make sure\n",
    "no overlapping patch is missed.\n",
    "This may well end up identifying a few patches that do not, in the end,\n",
    "overlap with the extent of the desired cutout -- and that's ok."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e15ec10-fc20-4e85-a3cb-c273d1fa0920",
   "metadata": {},
   "outputs": [],
   "source": [
    "sra1 = ra1 + box_side_deg\n",
    "sra2 = ra2 - box_side_deg\n",
    "sdec1 = dec1 - box_side_deg\n",
    "sdec2 = dec2 + box_side_deg\n",
    "\n",
    "scorners = np.asarray([[sra1, sdec1], [sra2, sdec1],\n",
    "                       [sra2, sdec2], [sra1, sdec2]],\n",
    "                      dtype='float')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c63856b-8d40-4a07-bc2b-e077f8e772b1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-04T21:16:07.067305Z",
     "iopub.status.busy": "2024-05-04T21:16:07.066559Z",
     "iopub.status.idle": "2024-05-04T21:16:07.070902Z",
     "shell.execute_reply": "2024-05-04T21:16:07.070415Z",
     "shell.execute_reply.started": "2024-05-04T21:16:07.067285Z"
    }
   },
   "source": [
    "Create the TAP query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "417e9968-15d1-4ea3-a336-4ab6317c8027",
   "metadata": {},
   "outputs": [],
   "source": [
    "string_polygon = \"\"\n",
    "for c, scorner in enumerate(scorners):\n",
    "    s1 = str(np.round(scorner[0], 4))\n",
    "    s2 = str(np.round(scorner[1], 4))\n",
    "    string_polygon += s1 + ', ' + s2\n",
    "    if c < 3:\n",
    "        string_polygon += ', '\n",
    "    del s1, s2\n",
    "\n",
    "query = \"SELECT lsst_patch,lsst_tract,s_dec,s_ra \"\\\n",
    "        \"FROM dp02_dc2_catalogs.CoaddPatches \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', s_ra, s_dec), \"\\\n",
    "        \"POLYGON('ICRS', \"+string_polygon+\"))=1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ea9cb99-03d8-4a3f-bc71-956e9b047dc3",
   "metadata": {},
   "source": [
    "Run the TAP query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e62d2bd-c496-4793-9ea0-ec1c9dfee7f0",
   "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": "markdown",
   "id": "c1db9b34-febc-4734-ad8d-c5521cd56dfd",
   "metadata": {},
   "source": [
    "If the returned job phase was \"error\", uncomment the following to review the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abba2e71-1a6b-484c-aae2-1c4bafb29bcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# job.raise_if_error()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d974a58-6878-4ee1-9cde-17629c11ef13",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-04T22:05:42.810839Z",
     "iopub.status.busy": "2024-05-04T22:05:42.810107Z",
     "iopub.status.idle": "2024-05-04T22:05:42.813866Z",
     "shell.execute_reply": "2024-05-04T22:05:42.813343Z",
     "shell.execute_reply.started": "2024-05-04T22:05:42.810818Z"
    }
   },
   "source": [
    "Retrieve the TAP results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "455e63aa-8aa1-4be4-b5a7-7d7aeccca5ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8afd1493-7625-4698-9e48-442920481213",
   "metadata": {},
   "source": [
    "Option to show the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d43c7da1-70e8-4581-9ad6-ffb291ff6c6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4508c32b-ed56-49b6-ace8-39f7d105ac5a",
   "metadata": {},
   "source": [
    "Hold the tract and patch values for images that \n",
    "will contribute, in addition to the anchor image,\n",
    "to the final custom cutout in `add_tracts_and_patches`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc5bfcc1-d755-448f-a247-d9e7854d1811",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = []\n",
    "for result in results:\n",
    "    tract = result['lsst_tract']\n",
    "    patch = result['lsst_patch']\n",
    "    if (tract != anchor_tract) | (patch != anchor_patch):\n",
    "        temp.append([tract, patch])\n",
    "add_tracts_and_patches = np.asarray(temp, dtype='int')\n",
    "del temp\n",
    "\n",
    "print('add these to the anchor patch: ')\n",
    "print(add_tracts_and_patches)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58ddf4b6-2e91-4c87-82fd-4d02af8db7a2",
   "metadata": {},
   "source": [
    "## 3. Create the custom cutout\n",
    "\n",
    "Define the filter of the desired custom cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "712f4871-69e3-46b0-b3fe-cb3f5c20d229",
   "metadata": {},
   "outputs": [],
   "source": [
    "use_filter = 'r'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de7b519-a918-4297-ab5a-19017ae18e51",
   "metadata": {},
   "source": [
    "### 3.1. Define and display the anchor image\n",
    "\n",
    "Define the `dataId` for the anchor image.\n",
    "\n",
    "Use the butler to retrieve the anchor image and its World Coordinate System (WCS)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6c45d64-ced7-4ee9-91a7-4aafc7ea07c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "anchor_dataId = {'tract': anchor_tract, 'patch': anchor_patch,\n",
    "                 'band': use_filter}\n",
    "anchor_image = butler.get('deepCoadd_calexp', dataId=anchor_dataId)\n",
    "anchor_wcs = butler.get('deepCoadd_calexp.wcs', dataId=anchor_dataId)\n",
    "print(anchor_wcs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c970cb2-6346-4c16-8499-d104cc063fac",
   "metadata": {},
   "source": [
    "Visually inspect the `anchor_image`. The below demonstration plots a circle centered at `my_ra` and `my_dec`, with diameter equal to the desired image box edge outlined above. The map indicates that the desired location is near the patch and tract edge and joining to other deepCoadds is required in order to display a contiguous image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fdb5624-712b-46c0-952a-d6118dbf23bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = box_side_deg/2. * u.deg\n",
    "coord = SkyCoord(ra=my_ra*u.degree, dec=my_dec*u.degree, frame='icrs')\n",
    "\n",
    "plt.subplot(projection=WCS(anchor_image.getWcs().getFitsMetadata()))\n",
    "extent = (anchor_image.getBBox().beginX, anchor_image.getBBox().endX,\n",
    "          anchor_image.getBBox().beginY, anchor_image.getBBox().endY)\n",
    "\n",
    "plt.imshow(anchor_image.image.array, vmin=.1, vmax=1, extent=extent,\n",
    "           origin='lower', cmap='gray')\n",
    "\n",
    "anchor_bbox = anchor_image.getBBox()\n",
    "\n",
    "aperture = SkyCircularAperture(coord, r=radius)\n",
    "pix_aperture = aperture.to_pixel(WCS(anchor_image.getWcs().getFitsMetadata()))\n",
    "pix_aperture.plot(color='r', lw=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a2bfcc9-448f-4f33-8a0b-9d4e3355e921",
   "metadata": {},
   "source": [
    "> Figure 2: The anchor image patch, with a partial red circle indicating that the size of the desired large cutout is bigger than a single `deepCoadd` patch."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9611d861-f18c-4a44-8be6-f11f93036332",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-04T22:11:37.852011Z",
     "iopub.status.busy": "2024-05-04T22:11:37.851584Z",
     "iopub.status.idle": "2024-05-04T22:11:37.875109Z",
     "shell.execute_reply": "2024-05-04T22:11:37.874640Z",
     "shell.execute_reply.started": "2024-05-04T22:11:37.851994Z"
    }
   },
   "source": [
    "### 3.2. Define the other images to combine\n",
    "\n",
    "Create the arrays of `dataId`s and images to be passed to `GetTemplateTask`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e2ab8a5-8d2e-4e71-9239-673bd9957e0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_dataIds = [anchor_dataId]\n",
    "all_images = [anchor_image]\n",
    "all_tracts = [anchor_tract]\n",
    "\n",
    "for i in range(len(add_tracts_and_patches)):\n",
    "    temp_tract = add_tracts_and_patches[i][0]\n",
    "    temp_patch = add_tracts_and_patches[i][1]\n",
    "    temp_dataId = {'tract': temp_tract, 'patch': temp_patch,\n",
    "                   'band': use_filter}\n",
    "    temp_image = butler.get('deepCoadd_calexp', dataId=temp_dataId)\n",
    "\n",
    "    all_images.append(temp_image)\n",
    "    all_dataIds.append(temp_dataId)\n",
    "    all_tracts.append(temp_tract)\n",
    "\n",
    "    del temp_tract, temp_patch, temp_dataId, temp_image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a08998f1-0cf7-46dd-8873-3ef2b83a03e7",
   "metadata": {},
   "source": [
    "Print the `dataId`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f71621c-c026-4b3e-8161-1eb20b4a5c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "for dataId in all_dataIds:\n",
    "    print(dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4ea673a-75e8-4322-aef9-0c1be00cdedc",
   "metadata": {},
   "source": [
    "### 3.3. Set up `getTemplateTask`\n",
    "\n",
    "Define `getTemplateTask` which will enable the various images to be warped and resampled into one big cutout. \n",
    "\n",
    "No configuration changes are needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c2cfe9d-5d8f-4e58-8dae-58b8120d7865",
   "metadata": {},
   "outputs": [],
   "source": [
    "getTemplateTask = GetTemplateTask()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33fe4b2f-943a-4bed-a189-948a980662fb",
   "metadata": {},
   "source": [
    "Build the WCS definition for the new cutout. There are two cases for how to proceed. \n",
    "\n",
    "#### 3.3.1. Within a single tract (limited use case)\n",
    "\n",
    "If the desired cutout crosses patch boundaries but is within a single tract, then the pixel grids are exactly compatible, and no reprojection is needed. Pixels from adjacent patches can be exactly stitched together at the boundaries, without having to re-define a WCS and reproject the pixels onto it.  Pixels can be taken from the relevant patches and stitched into a single array.  If not reprojecting, be careful to preserve the CRPIX/CRVAL values of the inputs and not attempt to recenter the WCS on the cutout, because it's still on a tangent plane relative to the tract center.\n",
    "\n",
    "This is _not_ the case that is demonstrated in this notebook, because the process for creating a big cutout that crosses tract boundaries is more general, and would work even if a tract boundary is not crossed.\n",
    "\n",
    "#### 3.3.2. Crossing tract boundaries (general use case)\n",
    "\n",
    "If the cutout crosses tract boundaries (which have independent projections and WCS), then it is necessary to reproject the input images to a new tangent plane centered at the target and with north up. This more general case is demonstrated below, and will work also for the simpler case of adjacent patches within one tract.\n",
    "\n",
    "The cell below defines the new WCS for the joined images across both patch and tract boundaries.\n",
    "\n",
    "The WCS consists of 4 components.\n",
    "\n",
    "1) The RA and Dec of the center of the image defined as `crval`.\n",
    "\n",
    "2) The center pixel (x, y) defined as `crpix`. <br>\n",
    "In Section 1, the cutout area of interest is defined in `box_side_deg` as 0.2 degrees on a side (12 arcminutes).\n",
    "The center is defined as `box_side_deg`, converted to arcsec and then divided in half.\n",
    "Dividing that by the `deepCoadd` pixel scale of 0.2 arcsec/pixel results in the center pixel, in units of pixels.\n",
    "\n",
    "3) The pixel scale, which is a 2x2 matrix defined as `cdMatrix`.<br>\n",
    "Note that the pixel scale of all deepCoadds are always the same, so the new WCS can simply use the one from `anchor_image`.\n",
    "\n",
    "5) The projection of the WCS `projection`, which is defined as `TAN` (referring to a new tangent projection centered at the `crval`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d648ade4-987d-4ec4-b9a5-e44bae5af41f",
   "metadata": {},
   "outputs": [],
   "source": [
    "crval = geom.SpherePoint(my_ra, my_dec, geom.degrees)\n",
    "\n",
    "crpix = box_side_deg * 3600 / 2. / 0.2\n",
    "crpix1 = crpix\n",
    "crpix2 = crpix\n",
    "crpix = geom.Point2D(crpix1, crpix2)\n",
    "\n",
    "x = anchor_image.getWcs().getFitsMetadata()\n",
    "cdMatrix = [x['CD1_1'], x['CD1_2']], [x['CD2_1'], x['CD2_2']]\n",
    "\n",
    "projection = 'TAN'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa831407-bf87-431d-abca-7cf38420b37f",
   "metadata": {},
   "source": [
    "Define the new WCS for the large cutout using the `makeSkyWcs` method based on the definition above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "299762b8-095b-4b0b-b36d-e3a91377d7da",
   "metadata": {},
   "outputs": [],
   "source": [
    "newWCS = makeSkyWcs(crpix, crval, cdMatrix, projection)\n",
    "print(newWCS)\n",
    "print(anchor_image.getWcs())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42cdee00-072e-4377-a060-24ce493eae3b",
   "metadata": {},
   "source": [
    "Create a bounding box for the new image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "838b76bb-5235-4caf-80a8-5651c8853d18",
   "metadata": {},
   "outputs": [],
   "source": [
    "corner_SE_deg = geom.SpherePoint(corners[0][0], corners[0][1], geom.degrees)\n",
    "corner_NW_deg = geom.SpherePoint(corners[2][0], corners[2][1], geom.degrees)\n",
    "corner_SE_pix = newWCS.skyToPixel(corner_SE_deg)\n",
    "corner_NW_pix = newWCS.skyToPixel(corner_NW_deg)\n",
    "\n",
    "xmin = int(np.floor(corner_SE_pix[0]))\n",
    "ymin = int(np.floor(corner_SE_pix[1]))\n",
    "xmax = int(np.floor(corner_NW_pix[0]))\n",
    "ymax = int(np.floor(corner_NW_pix[1]))\n",
    "\n",
    "newimg_shape = Image(geom.Box2I(minimum=geom.Point2I(x=xmin, y=ymin),\n",
    "                                maximum=geom.Point2I(x=xmax, y=ymax)),\n",
    "                     dtype=np.float32)\n",
    "newimage_bbox = newimg_shape.getBBox()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aab6d458-5922-4df8-94df-15fcc84ef4b2",
   "metadata": {},
   "source": [
    "### 3.4. Create the big cutout image\n",
    "\n",
    "Create the custom cutout as `newimage`.\n",
    "\n",
    ">**Notice:** Some input patch images are skipped for not overlapping with the defined bounding box after all, as was expected (Section 2.3). The argument `physical_filter` is required to be passed, even though use of r-band was already specified above.\n",
    "\n",
    "`getTemplateTask` expects the input lists of images and dataIds to be in dicts indexed on tract id.\n",
    "Use `collections.defaultdict(list)` to contain these, looping over unique tract values to create the per-tract lists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74dc1cd6-56df-4837-8fbf-13537fdfebff",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_dataIds_dict = collections.defaultdict(list)\n",
    "all_images_dict = collections.defaultdict(list)\n",
    "\n",
    "for tract in set(all_tracts):\n",
    "    dataId_list = []\n",
    "    image_list = []\n",
    "    for i in range(len(all_dataIds)):\n",
    "        if all_dataIds[i]['tract'] == tract:\n",
    "            dataId_list.append(all_dataIds[i])\n",
    "            image_list.append(all_images[i])\n",
    "    all_dataIds_dict[tract] = dataId_list\n",
    "    all_images_dict[tract] = image_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98a59774-18d1-4b23-9f29-f7dd3908d5c8",
   "metadata": {},
   "source": [
    "Print the dict of dataIds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2f00f24-549f-40f7-a453-f371854ef5db",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('dataIds', all_dataIds_dict)\n",
    "print('dict of images', all_images_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b5b5f91-d202-4b14-a562-275bb9c9ce76",
   "metadata": {},
   "outputs": [],
   "source": [
    "newimage = getTemplateTask.run(all_images_dict, newimage_bbox,\n",
    "                               newWCS, all_dataIds_dict,\n",
    "                               physical_filter='r_sim_1.4')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fd8e6cc-8289-42bc-a7b7-967b96ef1d9e",
   "metadata": {},
   "source": [
    "## 4. Visualize the custom cutout\n",
    "\n",
    "Extract metadata for the custom cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b1b7823-e3f2-480f-afb7-a189bbd60087",
   "metadata": {},
   "outputs": [],
   "source": [
    "newimage_wcs = newimage.template.getWcs()\n",
    "newimage_wcs_fmd = newimage.template.getWcs().getFitsMetadata()\n",
    "\n",
    "newimage_extent = (newimage_bbox.beginX, newimage_bbox.endX,\n",
    "                   newimage_bbox.beginY, newimage_bbox.endY)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27d8fcfb-a843-48ef-943b-aeb2a2ce542a",
   "metadata": {},
   "source": [
    "### 4.1. Display with matplotlib\n",
    "\n",
    "Compare the red circle in the larger cutout to that from the `anchor_image` in Section 3.1.\n",
    "The full desired area is now accessible in one single contiguous cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b1a16cb-19f3-46a6-816e-5436c06be048",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 8))\n",
    "plt.subplot(projection=WCS(newimage_wcs_fmd))\n",
    "plt.imshow(newimage.template.image.array, vmin=0.01, vmax=0.5,\n",
    "           extent=newimage_extent, origin='lower', cmap='gray')\n",
    "\n",
    "aperture = SkyCircularAperture(coord, r=radius)\n",
    "pix_aperture = aperture.to_pixel(WCS(newimage_wcs_fmd))\n",
    "pix_aperture.plot(color='r', lw=3)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64f75715-9f34-4f50-9fa1-1261df8f0fb0",
   "metadata": {},
   "source": [
    "> Figure 3: The new image (large cutout) with a red circle marking the user-defined size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "554f4f2f-a81f-4955-8c75-c3baeaf9f9b8",
   "metadata": {},
   "source": [
    "#### 4.1.1. Draw input image boundaries\n",
    "\n",
    "It can be useful to visualize and identify which region of the big joined image originated from which patch and tract.\n",
    "\n",
    "Below is an example demonstrating how to identify the origin image for pixels in the big image, based on the input list of `dataIds` that were identified in Section 2. \n",
    "\n",
    "The code runs through all of the `dataIds` of the neighboring patchs and tracts, identifies the pixel values of the image corners in the new image, and plots the rectangular regions defining the original images (each as a different color)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e4c1e38-2375-4a70-96d2-cbcaa299a445",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 8))\n",
    "plt.subplot(projection=WCS(newimage_wcs_fmd))\n",
    "plt.imshow(newimage.template.image.array, vmin=0.01, vmax=0.5,\n",
    "           extent=newimage_extent, origin='lower', cmap='gray')\n",
    "\n",
    "colors = ['red', 'orange', 'yellow', 'green', 'blue', 'purple']\n",
    "i = 0\n",
    "\n",
    "for dataId in all_dataIds:\n",
    "    print(dataId, colors[i])\n",
    "    dataId_image = butler.get('deepCoadd_calexp', dataId=dataId)\n",
    "    dataId_wcs = butler.get('deepCoadd_calexp.wcs', dataId=dataId)\n",
    "    dataId_bbox = dataId_image.getBBox()\n",
    "\n",
    "    corners2 = [dataId_wcs.pixelToSky(dataId_bbox.beginX, dataId_bbox.beginY),\n",
    "                dataId_wcs.pixelToSky(dataId_bbox.beginX, dataId_bbox.endY),\n",
    "                dataId_wcs.pixelToSky(dataId_bbox.endX, dataId_bbox.endY),\n",
    "                dataId_wcs.pixelToSky(dataId_bbox.endX, dataId_bbox.beginY)]\n",
    "\n",
    "    corners_dataId_sky = [geom.SpherePoint(corners2[0].getRa().asDegrees(),\n",
    "                                           corners2[0].getDec().asDegrees(),\n",
    "                                           geom.degrees),\n",
    "                          geom.SpherePoint(corners2[1].getRa().asDegrees(),\n",
    "                                           corners2[1].getDec().asDegrees(),\n",
    "                                           geom.degrees),\n",
    "                          geom.SpherePoint(corners2[2].getRa().asDegrees(),\n",
    "                                           corners2[2].getDec().asDegrees(),\n",
    "                                           geom.degrees),\n",
    "                          geom.SpherePoint(corners2[3].getRa().asDegrees(),\n",
    "                                           corners2[3].getDec().asDegrees(),\n",
    "                                           geom.degrees)]\n",
    "    corners_dataId_pix = newWCS.skyToPixel(corners_dataId_sky)\n",
    "\n",
    "    coord_box = [[corners_dataId_pix[0].getX(), corners_dataId_pix[0].getY()],\n",
    "                 [corners_dataId_pix[1].getX(), corners_dataId_pix[1].getY()],\n",
    "                 [corners_dataId_pix[2].getX(), corners_dataId_pix[2].getY()],\n",
    "                 [corners_dataId_pix[3].getX(), corners_dataId_pix[3].getY()]]\n",
    "\n",
    "    # store the pixel values of the NE corner;\n",
    "    # to be used to study PSF variation in Section 5\n",
    "    if i == 0:\n",
    "        NE_anchor_pix = np.floor(coord_box[2])\n",
    "\n",
    "    coord_box.append(coord_box[0])  # repeat the first to create a closed loop\n",
    "    xs, ys = zip(*coord_box)  # create lists of x and y values\n",
    "\n",
    "    plt.plot(xs, ys, color=colors[i])\n",
    "    i = i+1\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da655462-10da-4d50-bb76-6e10b1091c7c",
   "metadata": {},
   "source": [
    "> Figure 4: The new image, with the bounding boxes of the six individual `deepCoadd` images which were identified as potential contributors to the big cutout shown in orange, purple, blue, red, green, and yellow."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c728097-f6a1-4c77-ad10-d7487b2136b0",
   "metadata": {},
   "source": [
    "### 4.2. Display with `afwDisplay`\n",
    "\n",
    "As an option, display the big coadd with `afwDisplay`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2b7ead9-a7bd-4cc2-9aca-0f441b9ff003",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 8))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(newimage.template.image)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b746719d-ad24-4285-985c-a64d80a105ca",
   "metadata": {},
   "source": [
    "> Figure 5: The new image, displayed in greyscale with `afwDisplay`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c749612b-8130-431a-9d7e-34951883f9b5",
   "metadata": {},
   "source": [
    "### 4.3. Display with Firefly\n",
    "\n",
    "Firefly allows for interactive image display and manipulation.\n",
    "See DP0.2 tutorial notebook 03b for a full demonstration of how to use Firefly within the Notebook Aspect.\n",
    "\n",
    "Reset the `afwDisplay` backend to be Firefly and start the display.\n",
    "A new tab will open, containing the Firefly interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7e26c4f-f426-4b5c-85fb-105f7090dede",
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('firefly')\n",
    "afw_display = afwDisplay.Display(frame=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af11f3d0-37fc-4041-986c-27d4ed726f11",
   "metadata": {},
   "source": [
    "Display the image in Firefly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15df02e5-3dd9-45fe-8c04-5d9e464a08ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display.mtv(newimage.template)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cf35ffd-8a3a-4941-9857-476b9796e342",
   "metadata": {},
   "source": [
    "The Firefly default is to visualize the mask plane with colors. \n",
    "Set the mask plane transparency to 100, fully transparent, to see only the pixel data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6af04896-ba12-42e3-9498-baf779d7361c",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display.setMaskTransparency(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "399da00b-9f89-4e83-b7cf-885d43df89e4",
   "metadata": {},
   "source": [
    "### 4.4. Option to save as a FITS file\n",
    "\n",
    "To save as a FITS file in the home directory, uncomment (remove the #)\n",
    "from all lines in the following cell and then execute it.\n",
    "\n",
    "> **Warning:** Images can take up a lot of disk space. Save and download with caution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8fa969f-70cf-4548-93ca-9db5bf9b217a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# username = os.environ.get('USER')\n",
    "# fnm = '/home/' + username + '/my_big_cutout.fits'\n",
    "# newimage.template.writeFits(fnm)\n",
    "# del username, fnm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e55954a6-27e7-432e-bb99-3453afe3c59b",
   "metadata": {},
   "source": [
    "## 5. Accessing the PSF metadata of the new cutout\n",
    "\n",
    "Where the tracts overlap, the pixel values are averaged to generate the resulting template image.\n",
    "\n",
    "The point-spread function (PSF) on the template (big cutout) is created by combining the `CoaddPsf` on each template image into a meta-CoaddPsf. Thus, it is possible to extract the PSF properties from any region of the new big cutout image, in the same way it is done for the `deepCoadd` images.\n",
    "\n",
    "Below demonstrates the way to access the PSF meta data from the template image.  \n",
    "\n",
    "More information on how to access PSF information can be found in DP0.2 tutorial notebooks 12a and 12b."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "323f1935-f20a-4de3-a151-9237ce79d152",
   "metadata": {},
   "source": [
    "### 5.1. Compare the PSF size\n",
    "\n",
    "First, check that the PSFs in the original image is the same as in that location of the template image. The first cell below will compare the PSFs at the location of the target, where the PSF is averaged across two patches. The shape of the PSF in the template is the same as the original to within 0.1-0.2%.\n",
    "\n",
    "\n",
    "Compare the size (the width; the standard deviation or $\\sigma$, `sigma`) of the PSF of the original (anchor) image\n",
    "to that of the new (template) image for a fixed sky location: the user-defined center of the new image (`center`).\n",
    "\n",
    "At this location, two patch images overlapped (as shown in Section 4.1.1.), so pixel data has been averaged, and some\n",
    "small variation is to be expected.\n",
    "\n",
    "Below it is shown that the PSF size in the anchor image is approximately the same in the new image, for the given location, \n",
    "within 0.2%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1a74162-14a8-4601-95e0-ca03adea39ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "center = geom.SpherePoint(my_ra, my_dec, geom.degrees)\n",
    "\n",
    "anchor_info = anchor_image.getInfo()\n",
    "center_pix_anchor = anchor_wcs.skyToPixel(center)\n",
    "newimage_info = newimage.template.getInfo()\n",
    "center_pix_newimage = newWCS.skyToPixel(center)\n",
    "point_tuple = (center_pix_anchor[0], center_pix_anchor[1])\n",
    "point_image = geom.Point2D(point_tuple)\n",
    "anchor_psf = anchor_image.getPsf()\n",
    "psf_shape = anchor_psf.computeShape(point_image)\n",
    "psf_sigma_old = psf_shape.getDeterminantRadius()\n",
    "print(\"anchor center sigma=\", psf_sigma_old)\n",
    "\n",
    "newimage_psf = newimage.template.getPsf()\n",
    "point_tuple = (center_pix_newimage[0], center_pix_newimage[1])\n",
    "point_image = geom.Point2D(point_tuple)\n",
    "psf_shape = newimage_psf.computeShape(point_image)\n",
    "psf_sigma_new = psf_shape.getDeterminantRadius()\n",
    "print(\"newimage center sigma=\", psf_sigma_new)\n",
    "\n",
    "percent_diff = np.abs(psf_sigma_old - psf_sigma_new)/psf_sigma_new * 100\n",
    "print(\"percent difference is \", percent_diff, \"%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "036a0ede-3a51-4403-ae72-d9af59c4b4f2",
   "metadata": {},
   "source": [
    "Next, check that the size of the PSF in the anchor image is the same as in the new image, for a location that is not averaged across multiple patches, but is warped to the new WCS.\n",
    "\n",
    "The cell below compares the size of the PSFs at the location of the southeast corner (`corner_SE_deg` as defined in Section 2), where the PSF is derived only from the single original patch.\n",
    "\n",
    "The size of the PSF in the new image is also the same as the anchor image to within $\\sim$0.2%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef768a21-71ce-44d3-aeab-12fe809f44b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "center_pix = anchor_wcs.skyToPixel(corner_SE_deg)\n",
    "point_tuple = (center_pix[0], center_pix[1])\n",
    "point_image = geom.Point2D(point_tuple)\n",
    "psf_shape = anchor_psf.computeShape(point_image)\n",
    "psf_sigma_old = psf_shape.getDeterminantRadius()\n",
    "print(\"anchor SE corner sigma=\", psf_sigma_old)\n",
    "\n",
    "center_pix = newWCS.skyToPixel(corner_SE_deg)\n",
    "point_tuple = (center_pix[0], center_pix[1])\n",
    "point_image = geom.Point2D(point_tuple)\n",
    "psf_shape = newimage_psf.computeShape(point_image)\n",
    "psf_sigma_new = psf_shape.getDeterminantRadius()\n",
    "print(\"newimage SE corner sigma=\", psf_sigma_new)\n",
    "\n",
    "percent_diff = np.abs(psf_sigma_old - psf_sigma_new)/psf_sigma_new * 100\n",
    "print(\"percent difference is \", percent_diff, \"%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "396f798a-b0e9-480f-9655-8521884fbc76",
   "metadata": {},
   "source": [
    "### 5.2. Spatial variation in PSF differences\n",
    "The following cell will more generally visualize the spatial variation of the percent difference of the PSF size using a heatmap, for just the part of anchor_image that overlaps the newimage.  To start, define the northwest and southwest corners of the overlapping region between `anchor_image` and `newImage` (`NE_anchor_pix` was stored in Section 4.1 for this purpose). Before display, the array has to be transposed using `.T` to match the orientation of the sky coordinates in the map above.\n",
    "\n",
    "The heat map of percent difference in PSF sizes across the map indicate that any spatial variation is incredibly small, < 0.3 %. Note that the percent difference between the PSF of `anchor_image` and `newImage` increases at x-pixels > 1600, which marks the boundary where a second deepCoadd also contributes to the new image pixels and PSF properties (yellow box in Section 4.1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5708010d-eac5-4952-8dda-d27187f2bf47",
   "metadata": {},
   "outputs": [],
   "source": [
    "NW_x = NE_anchor_pix[0]\n",
    "NW_y = NE_anchor_pix[1]\n",
    "\n",
    "SE_x = 0\n",
    "SE_y = 0\n",
    "\n",
    "x_arr = np.arange(SE_x, NW_x, 300)\n",
    "y_arr = np.arange(SE_y, NW_y, 300)\n",
    "\n",
    "PSF_diff = np.zeros((len(x_arr), len(y_arr)))\n",
    "\n",
    "for i in range(len(x_arr)):\n",
    "    for j in range(len(y_arr)):\n",
    "        point_new = geom.Point2D(x_arr[i], y_arr[j])\n",
    "        sigma_new = newimage_psf.computeShape(point_new).getDeterminantRadius()\n",
    "        coords4anchor = newWCS.pixelToSky(x_arr[i], y_arr[j])\n",
    "        point_old = geom.Point2D(anchor_wcs.skyToPixel(coords4anchor))\n",
    "        sigma_old = anchor_psf.computeShape(point_old).getDeterminantRadius()\n",
    "        PSF_diff[i][j] = np.abs(sigma_old - sigma_new) / sigma_new * 100\n",
    "\n",
    "plt.imshow(PSF_diff.T, cmap='rainbow',\n",
    "           extent=[np.min(x_arr), np.max(x_arr),\n",
    "                   np.min(y_arr), np.max(y_arr)])\n",
    "\n",
    "plt.colorbar(label='Percent difference in PSF size')\n",
    "plt.xlabel('new image pixels')\n",
    "plt.ylabel('new image pixels')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98d531b2-68a3-4b4d-b2a1-68980f2b8748",
   "metadata": {},
   "source": [
    "> Figure 6: A heatmap showing the difference in PSF size between the anchor image and the new image, only for the area of the new image to which the anchor image contributed. The colorbar ranges from 0.12% (purple) to 0.26% (red). The difference in PSF size varies smoothly in the left half, where the new image is composed just from the anchor image. However, the difference in PSF size varies randomly in the right half, where the new image is composed of pixel data from both the anchor image and an adjacent image."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65cab174-33d0-43c5-b8cc-39232421f363",
   "metadata": {},
   "source": [
    "## Appendix: Alternative solution using the getOverlappingExposures method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c598fecb-93e5-4a8d-af58-153376e3b797",
   "metadata": {},
   "source": [
    "The identification of the patch/tract images that overlap the big cutout area of interest can also be performed using the `getOverlappingExposures` method, without needing to build the new image using `getTemplateTask.run`. This method will take as input the LSST images that may overlap the desired area in the format of a `DeferredDatasetHandle`, and return a dictionary of lists containing those only of images that overlap that area. While the exclusion of non-overlapping images is done automatically by `getTemplateTask.run`, the example below shows how to use `getOverlappingExposures` to identify the overlapping images without having to build the big cutout with the run method (which will only print the overlapping images, but not store them). \n",
    "\n",
    "Further information on `datasetRef` types can be found in DP0.2 Tutorial Notebook 04b.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4edf8871-8b97-4bb9-a73a-2722a06ae48f",
   "metadata": {},
   "source": [
    "The `getOverlappingExposures` method of `getTemplateTask` is documented <a href=\"https://pipelines.lsst.io/v/daily/py-api/lsst.ip.diffim.GetTemplateTask.html#lsst.ip.diffim.GetTemplateTask.getOverlappingExposures\">here</a> and requires as input a dictionary that includes at least 4 components: a list of lsst `ExposureF` objects in `DeferredDatasetHandle` format; the template bounding box for the desired big coadd (as defined in section 3); the input geometry (SkyMap) for the template exposures, and the template WCS onto which the coadds will be resampled (also as defined in section 3). We provide a simple example below that does the same thing as Section 3 as a demonstration of using the output from `getOverlappingExposures`.\n",
    "\n",
    "The cell below demonstrates how to get a list of DeferredDatasetHandle datarefs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "051fb7d4-4772-468c-a45c-1dcec6c507ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "dI = all_dataIds[0]\n",
    "dataset_refs = list(butler.query_datasets('deepCoadd', data_id=dI))\n",
    "dataset_refs = set(dataset_refs)\n",
    "\n",
    "for dI in all_dataIds[1:]:\n",
    "    dr = list(butler.query_datasets('deepCoadd', data_id=dI))\n",
    "    dataset_refs = dataset_refs.union(set(dr))\n",
    "\n",
    "coadd_deferred_dataset_handle = [butler.getDeferred(dr) for dr in dataset_refs]\n",
    "sky_map = butler.get(\"skyMap\")\n",
    "inputs = {\"coaddExposures\": coadd_deferred_dataset_handle,\n",
    "          \"bbox\": newimage_bbox,\n",
    "          \"skyMap\": sky_map,\n",
    "          \"wcs\": newWCS, }\n",
    "\n",
    "print(\"Collecting the patches/tracts that overlap the desired bounding box\")\n",
    "results = getTemplateTask.getOverlappingExposures(inputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "776c80c8-2978-42ef-a3b9-6849a19df1a7",
   "metadata": {},
   "source": [
    "Below, access the metadata about the overlapping exposures (noting that the output is indexed by tract):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aa9d195-e2ff-424d-aa8e-cc965511252b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('dictionary of dataIds for overlapping exposures:', results.dataIds)\n",
    "print('dictionary of overlapping exposures:', results.coaddExposures)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bcb657d-e8a1-433f-ac6f-8cb06804d10d",
   "metadata": {},
   "source": [
    "Optionally, generate the big cutout image based on these overalpping exposures (as in section 3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6eef32e6-34fe-49cd-9d5a-8b0d45fc4d55",
   "metadata": {},
   "outputs": [],
   "source": [
    "newimage_appendix = getTemplateTask.run(results.coaddExposures, newimage_bbox,\n",
    "                                        newWCS, results.dataIds,\n",
    "                                        physical_filter='r_sim_1.4')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "391043c5-2dc8-4456-96d6-e47f7b687d8c",
   "metadata": {},
   "source": [
    "Optionally, plot the newly generated big cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "028f2e64-d9d5-4cc0-a3ad-b440f600473c",
   "metadata": {},
   "outputs": [],
   "source": [
    "newimanewimage_wcs = newimage_appendix.template.getWcs()\n",
    "newimage_wcs_fmd = newimage_appendix.template.getWcs().getFitsMetadata()\n",
    "\n",
    "newimage_extent = (newimage_bbox.beginX, newimage_bbox.endX,\n",
    "                   newimage_bbox.beginY, newimage_bbox.endY)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "plt.subplot(projection=WCS(newimage_wcs_fmd))\n",
    "plt.imshow(newimage_appendix.template.image.array, vmin=0.01, vmax=0.5,\n",
    "           extent=newimage_extent, origin='lower', cmap='gray')\n",
    "\n",
    "aperture = SkyCircularAperture(coord, r=radius)\n",
    "pix_aperture = aperture.to_pixel(WCS(newimage_wcs_fmd))\n",
    "pix_aperture.plot(color='r', lw=3)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}