DP02_10_Deblender_Data_Products.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# Introduction to the LSST deblender data products\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br>\n",
    "Contact author(s): Christina Williams <br>\n",
    "Last verified to run: 2024-10-29 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_42 <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 tutorial introduces the data products created by the LSST multiband deblender (Scarlet; <a href=\"https://ui.adsabs.harvard.edu/abs/2018A%26C....24..129M/abstract\">Melchior et al. 2018</a>) and illustrates some simple use cases for them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use of the catalog data products related to deblending objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** images (deepCoadd, skyMap, MaskedImage) and catalogs (objectTable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** lsst.afw.image, lsst.afw.detection, lsst.rsp, lsst.daf.butler, lsst.geom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:** Developed by Christina Willliams in collaboration with Melissa Graham and the Rubin Community Science Team for DP0.2.\n",
    "Influenced by Stack Club notebooks demonstrating how to run the LSST deblender (Scarlet; Melchior et al. 2018) by Fred Moolekamp and Alex Drlica-Wagner. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28e91cbf-ab7f-4e26-9276-b00299d6065e",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. Questions are welcome as new topics in the <a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> of the Rubin Community Forum. Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook is an introduction to the data products that are associated with the deblending process.\n",
    "\n",
    "In order to build its source catalog and measure photometry of objects in images, the LSST Science Pipeline uses the multi-wavelength deblending algorithm Scarlet (<a href=\"https://ui.adsabs.harvard.edu/abs/2018A%26C....24..129M/abstract\">Melchior et al. 2018</a>) to identify independent detected components within source blends and separate the flux from independent objects.\n",
    "\n",
    "A \"**parent**\" source is identified as a region of connected pixels above a given signal-to-noise ratio (S/N) in the `deepCoadd` image.\n",
    "These regions are called footprints.\n",
    "Each footprint may have one or more peaks, and it is these peaks that the deblender will use to infer the number and positions of objects blended in each footprint.\n",
    "The deblended sub-peaks are referred to as \"**children**\" of the parent in the properties saved by the deblender in the `Object` table.\n",
    "\n",
    "The Scarlet deblender separates the flux among components by modeling the flux and shape of each source making up the blend.\n",
    "\n",
    "This tutorial will identify a heavily blended parent object with many children and explore how to interpret the deblender-related measurements in the `Object` table.\n",
    "It will also demonstrate how to use deblender-related columns in the `Object` table to identify unique deblended child sources, and\n",
    "provide an example of how to access metadata that is available for each deblended source (e.g. the deblended footprints and pixel-weight maps).\n",
    "\n",
    "#### Additional resources\n",
    "\n",
    "An introduction to all of the available deblending flags, and the processes that produce the deblender data products discussed in this tutorial\n",
    "can be found in this <a href=\"https://pipelines.lsst.io/modules/lsst.pipe.tasks/deblending-flags-overview.html\">overview of the deblending flags</a>.\n",
    "\n",
    "A general discussion of blending impacts to LSST can be found in this article titled \"<a href=\"https://lss.fnal.gov/archive/2021/pub/fermilab-pub-21-598-ppd.pdf\">The challenge of blending in large sky surveys</a>\" by Melchior et al. 2021. \n",
    "\n",
    "Some more in-depth presentations on the deblending implementation in the LSST Science Pipelines are available from recorded talks during the Rubin Project and Community Workshop 2022's session on \"<a href=\"https://project.lsst.org/meetings/rubin2022/agenda/deblending-plans-and-challenges\">Deblending Plans and Challenges</a>\". "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1. Package Imports\n",
    "\n",
    "The [matplotlib](https://matplotlib.org/) (and especially sublibrary `matplotlib.pyplot`), [numpy](http://www.numpy.org/), and [astropy](http://www.astropy.org/) libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis.\n",
    "\n",
    "The `lsst.rsp` package provides access to the Table Access Protocol (TAP) service for queries to the DP0 catalogs.\n",
    "\n",
    "The `lsst.afw.image` and `lsst.afw.detection` packages provide access to some of the extended products created by the deblender pipeline. \n",
    "\n",
    "The `lsst.afw.display` library provides access to image visualization routines and the `lsst.daf.butler` library is used to access data products via the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d185358a-4c15-4b04-8e3d-13083191ee72",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import lsst.daf.butler as Butler\n",
    "from lsst.rsp import get_tap_service\n",
    "import lsst.geom as geom\n",
    "import lsst.afw.display as afwDisplay\n",
    "from lsst.afw.image import MaskedImage"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters\n",
    "\n",
    "The following two functions enable the visualization of deblended products:\n",
    " 1. `cutout_coadd` will make an image cutout (see DP0.2 tutorial notebook 03a on image display).\n",
    " 2. `heavyFootprint2Image` will help us grab the \"footprint\" of deblended objects for further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80eef068-b4db-46de-80fe-539df61c642a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def cutout_coadd(butler, ra, dec, band='r', datasetType='deepCoadd',\n",
    "                 skymap=None, cutoutSideLength=51, **kwargs):\n",
    "    \"\"\"\n",
    "    Produce a cutout from a coadd at the given ra, dec position.\n",
    "\n",
    "    Adapted from DC2 tutorial notebook by Michael Wood-Vasey.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Helper object providing access to a data repository\n",
    "    ra: float\n",
    "        Right ascension of the center of the cutout, in degrees\n",
    "    dec: float\n",
    "        Declination of the center of the cutout, in degrees\n",
    "    band: string\n",
    "        Filter of the image to load\n",
    "    datasetType: string ['deepCoadd']\n",
    "        Which type of coadd to load.  Doesn't support 'calexp'\n",
    "    skymap: lsst.afw.skyMap.SkyMap [optional]\n",
    "        Pass in to avoid the Butler read.  Useful if you have lots of them.\n",
    "    cutoutSideLength: float [optional]\n",
    "        Size of the cutout region in pixels.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \"\"\"\n",
    "    radec = geom.SpherePoint(ra, dec, geom.degrees)\n",
    "    cutoutSize = geom.ExtentI(cutoutSideLength, cutoutSideLength)\n",
    "\n",
    "    if skymap is None:\n",
    "        skymap = butler.get(\"skyMap\")\n",
    "\n",
    "    tractInfo = skymap.findTract(radec)\n",
    "    patchInfo = tractInfo.findPatch(radec)\n",
    "    xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))\n",
    "    bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)\n",
    "    print(\"bbox = \", bbox, \"xy = \", xy, \"cutoutSize = \", cutoutSize)\n",
    "    patch = tractInfo.getSequentialPatchIndex(patchInfo)\n",
    "\n",
    "    coaddId = {'tract': tractInfo.getId(), 'patch': patch, 'band': band}\n",
    "    parameters = {'bbox': bbox}\n",
    "\n",
    "    cutout_image = butler.get(datasetType, parameters=parameters,\n",
    "                              dataId=coaddId)\n",
    "\n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba97da22-c4b8-480b-9792-39f9800c7b35",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def heavyFootprint2Image(heavy, fill=np.nan, bbox=None, imageType=MaskedImage):\n",
    "    \"\"\"\n",
    "    Create an image of a HeavyFootprint\n",
    "    Written by Fred Moolekamp\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    heavy : `HeavyFootprint`\n",
    "        The HeavyFootprint to insert into the image\n",
    "    fill: number\n",
    "        Number to fill the pixels in the image that are not\n",
    "        contained in `heavy`.\n",
    "    bbox : `Box2I`\n",
    "        Bounding box of the output image.\n",
    "    imageType : `type`\n",
    "        This should be either a `MaskedImage` or `Image` and describes\n",
    "        the type of the output image.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    image : `lsst.afw.image.MaskedImage` or `lsst.afw.image.Image`\n",
    "        An image defined by `bbox` and padded with `fill` that\n",
    "        contains the projected flux in `heavy`.\n",
    "    \"\"\"\n",
    "    if bbox is None:\n",
    "        bbox = heavy.getBBox()\n",
    "    image = imageType(bbox, dtype=heavy.getImageArray().dtype)\n",
    "    image.set(fill)\n",
    "    heavy.insert(image)\n",
    "    return image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a894b55a-354b-4108-b6da-097ccd3a5970",
   "metadata": {},
   "source": [
    "Define parameters for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54e8800e-26f3-4d70-9306-ef3a5fd75aa9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "plt.style.use('tableau-colorblind10')\n",
    "\n",
    "params = {'axes.labelsize': 28,\n",
    "          'font.size': 24,\n",
    "          'legend.fontsize': 14,\n",
    "          'xtick.major.width': 3,\n",
    "          'xtick.minor.width': 2,\n",
    "          'xtick.major.size': 12,\n",
    "          'xtick.minor.size': 6,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 3,\n",
    "          'axes.linewidth': 3,\n",
    "          'axes.labelweight': 3,\n",
    "          'axes.titleweight': 3,\n",
    "          'ytick.major.width': 3,\n",
    "          'ytick.minor.width': 2,\n",
    "          'ytick.major.size': 12,\n",
    "          'ytick.minor.size': 6,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [8, 8],\n",
    "          'figure.facecolor': 'White'}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd5b5811-70fd-4b9a-8d60-5144943ebc9f",
   "metadata": {},
   "source": [
    "Start the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ddfb55f-7810-4238-a287-877b79a84d78",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")\n",
    "assert service is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d59e33-12cd-4eb0-9cd9-02dfda57d37d",
   "metadata": {},
   "source": [
    "Instantiate the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17f8537b-387c-4974-b5e4-148ebfb9cfbc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "butler = Butler.Butler('dp02', collections='2.2i/runs/DP0.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec51ac0b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. The deblender data products\n",
    "\n",
    "The object and source catalogs contain a number of columns that can be used to identify and characterize blended sources (parents) that have been deblended into child sources, both of which are ultimately stored in the catalogs along with isolated single sources that do not require deblending.\n",
    "These columns help the user identify a unique sample and avoid returning duplicate objects in queries.\n",
    "\n",
    "In this tutorial we will focus on deblending data products for `Objects` detected and measured in the `deepCoadd` images.\n",
    "\n",
    "The following describes boolean flags that are set by the deblender, and some measurements or properties of deblended sources that can be used to characterize (de)blended sources.\n",
    "\n",
    "> **Note:** This tutorial uses catalog data that are available via TAP **and/or** the butler.\n",
    "In the DP0-era, some deblending data products are available *only* via butler-based catalogs, and cannot be accessed by TAP and are not listed in the <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html\">DP0.2 Schema Browser</a>.\n",
    "\n",
    "**Available in the TAP `Object` table and in the butler-based `objectTable`**:\n",
    "\n",
    "`detect_isPrimary` (boolean): True if source has no children, is in inner region of a coadd patch, and is in the inner region of a coadd tract (i.e., is a unique object in the catalog).  \n",
    "Note that an object where `detect_isPrimary`  = True can still be a deblended child, but itself was not a parent that was separated by the deblender. \n",
    "\n",
    "> **Warning: to avoid returning duplicates use `detect_isPrimary` = True for all catalog queries.**\n",
    "(Duplicates can be caused by including both the parent blend and the nearest child source, or by including the overlapping area outside the \"inner region\" of a patch or tract.\n",
    "\n",
    "`detect_fromBlend` (boolean): If True, this source is deblended from a parent with more than one child.\n",
    "\n",
    "`detect_isIsolated` (boolean): If True, this source was not part of a blend (i.e., the source has a single peak, `deblend_nPeaks` = 1, or its parent only had a single peak, `deblend_parentNPeaks` = 1). \n",
    "\n",
    "`detect_isDeblendedSource` (boolean): True if source has no children, is in the inner regions of its coadd patch and tract (will be unique) and is either an unblended isolated source or a deblended child from a parent.\n",
    "\n",
    "`detect_isDeblendedModelSource` (boolean): True if the source  has no children and is a deblended child (but this flag is false for unblended isolated sources); similar to `detect_isDeblendedSource`.\n",
    " \n",
    "`deblend_skipped` (boolean): If True, the deblender skipped this source and did not attempt to deblend.\n",
    "\n",
    "`parentObjectId` (long): Unique ID of parent source. \n",
    "\n",
    "`deblend_nChild` (int): Number of children this object has (defaults to 0).\n",
    "Differs from `deblend_nPeaks` in that `deblend_nChild` excludes both isolated sources and child peaks that were culled during deblending.\n",
    "\n",
    "`footprintArea` (int): Number of pixels in the source's detection \"footprint\".\n",
    "The footprint itself is a boolean mask of pixels in the input image that were detected as part of the parent blend.\n",
    "This corresponds to the footprint for the \"reference band\".\n",
    "\n",
    "`refBand` (char): The filter that the deblender identified the \"best\" location for the centroid of the peak flux based on the multi-wavelength deblending routine (referred to as the \"reference filter\").\n",
    "\n",
    "`<f>_blendedness` (dbl): For each filter (f = u, g, r, i, z, or y), a measure of how much the flux is affected by neighbors, (1 - child\\_flux/parent\\_flux).\n",
    "Operates on the absolute value of the pixels to try to obtain a de-noised value.\n",
    "See section 4.9.11 of <a href=\"https://ui.adsabs.harvard.edu/abs/2018PASJ...70S...5B/abstract\">Bosch et al. 2018, PASJ, 70, S5</a> for details. \n",
    "\n",
    "**Available only in the butler-based `objectTable`**:\n",
    "\n",
    "`deblend_nPeaks` (int): Number of peaks contained in the parent blend's footprint (is 1 for isolated sources).\n",
    "\n",
    "`deblend_parentNPeaks` (int): Number of peaks of the deblended child's parent."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d8dd998-f58f-4cad-8747-4764768b82a3",
   "metadata": {},
   "source": [
    "### 2.1. Find parents with many children\n",
    "\n",
    "\n",
    "Use the TAP service to find a bunch of parent blended objects.\n",
    "\n",
    "As with all explorative TAP queries searching for examples, it is essential to first start with a small area while developing queries (to confirm the query is good), and then iteratively increase the area as needed (because the full query could take hours without any spatial constraint on the search). This is because the tables are indexed by coordinate (RA, Dec) and queries with constraints only on measurement columns\n",
    "(such as deblending-related parameters) are inefficient.\n",
    "In other words, Qserv (the backend hosting the TAP-accessible tables) is \"spatially sharded\".\n",
    "For more information about efficient TAP queries, see tutorial 02.\n",
    "\n",
    "Search in the center of the DC2 simulation for wide-area (`footprintArea` > 20000, or about 30\"x30\") parent objects with many children (`deblend_nChild` > 20) that are not point-source like (`refExtendedness` = 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7e280d4-1c81-4396-884a-e5a2cef3f105",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ra = 55.064\n",
    "dec = -29.783\n",
    "query = \"SELECT objectId, coord_ra, coord_dec, x, y, tract, patch, \" + \\\n",
    "        \"refExtendedness, detect_isPrimary, detect_fromBlend, \" + \\\n",
    "        \"deblend_nChild, footprintArea \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec ), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + str(ra) + \", \" + str(dec) + \", 0.2)) = 1 \" + \\\n",
    "        \"AND deblend_nChild > 20 AND footprintArea > 20000 \" + \\\n",
    "        \"AND refExtendedness = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eba8a24-a92f-4bfa-aa42-90a1135d597d",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85742354-e885-4176-957b-ca3c28562fde",
   "metadata": {},
   "outputs": [],
   "source": [
    "parents = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e31ffb14-3d91-4f39-8c9c-01a815589b50",
   "metadata": {},
   "outputs": [],
   "source": [
    "parents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88e3d39c-779e-461f-837e-23dbbf312033",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-02-28T01:40:02.089591Z",
     "iopub.status.busy": "2023-02-28T01:40:02.089206Z",
     "iopub.status.idle": "2023-02-28T01:40:02.095587Z",
     "shell.execute_reply": "2023-02-28T01:40:02.094787Z",
     "shell.execute_reply.started": "2023-02-28T01:40:02.089563Z"
    }
   },
   "source": [
    "Use the `cutout_coadd` function to obtain a cutout image of one of these parents (use parent index `pi=1`).\n",
    "Make the side length 3x the square root of the `footprintArea`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02bd02b8-cb35-4735-b5a2-a38b3501494c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pi = 1\n",
    "cutout = cutout_coadd(butler, parents['coord_ra'][pi],\n",
    "                      parents['coord_dec'][pi],\n",
    "                      band='i', datasetType='deepCoadd',\n",
    "                      cutoutSideLength=3.0\n",
    "                      * np.sqrt(float(parents['footprintArea'][pi])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cf08c3a-7244-4f43-a25c-a446f615bf36",
   "metadata": {},
   "source": [
    "Use TAP to retrieve all children for the selected parent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72a062c5-b808-4dee-ac74-504c90a4bf62",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, coord_ra, coord_dec, x, y \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec ), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + str(ra) + \", \" + str(dec) + \", 0.3)) = 1 \" + \\\n",
    "        \"AND parentObjectId = \" + str(parents['objectId'][pi]) + \" \" + \\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da2109a6-8a3c-42b0-9f9e-1e9d5cfbc2bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6047224-1d9c-4765-88ed-13fba3d37dc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "children = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d5198ba-d40f-44e8-a561-529101cae2c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# children"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a7a40da-dcae-45f1-b1c8-19c8544bdb11",
   "metadata": {},
   "source": [
    "Display the `deepCoadd` cutout made for this parent, with a blue cross for the parent and red circles for its children."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bae3bc26-68f3-41a4-a7e4-8624667b867d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout.image)\n",
    "with display.Buffering():\n",
    "    display.dot('+', parents['x'][pi], parents['y'][pi], ctype=afwDisplay.BLUE)\n",
    "    for ci in range(len(children)):\n",
    "        display.dot('o', children['x'][ci], children['y'][ci],\n",
    "                    size=10, ctype=afwDisplay.RED)\n",
    "plt.gca().axis('off')\n",
    "display.show_colorbar(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b019b8-2f36-49e3-b428-e48d5e5bd3f3",
   "metadata": {},
   "source": [
    "> Figure 1: The cutout image displayed in greyscale, with a blue + marking the parent object and red circles marking the deblended children. In this case, the blended object looks like a rich galaxy cluster.\n",
    "\n",
    "> **Note:** Performance of the deblender is still being optimized, and results like this within DP0 are not necessarily representative of future results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60fffdb0-8a69-47f5-85a6-d41c2095d09e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del parents, children, cutout"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4510ebd7-0c2e-456a-81d9-49fc28db381b",
   "metadata": {},
   "source": [
    "### 2.2. Find a parent at a given location\n",
    "\n",
    "In the case where the user wants to explore the (de)blended objects at a specific location, the process is similar to the one above except here, no constraints are made _except_ spatial constraints.\n",
    "\n",
    "In this case our search area is limited to a circle of radius 0.003 degrees, or about 10 arcseconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6430d99-187f-48fa-a9f1-af068e52c6d3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ra = 50.1089143\n",
    "dec = -44.4812763\n",
    "query = \"SELECT objectId, coord_ra, coord_dec, x, y, tract, patch, \" + \\\n",
    "        \"detect_isPrimary, deblend_nChild, deblend_skipped, \" + \\\n",
    "        \"detect_fromBlend, detect_isDeblendedModelSource, \" +\\\n",
    "        \"detect_isDeblendedSource, detect_isIsolated, parentObjectId, \" + \\\n",
    "        \"footprintArea, z_blendedness \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec ), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + str(ra) + \", \" + str(dec) + \", 0.003)) = 1 \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8777edb-0c9c-480c-8de3-d0463a5aaf86",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2150eb8-7882-44f5-8b93-304052da8420",
   "metadata": {},
   "outputs": [],
   "source": [
    "table = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61642f4a-4135-450b-baa2-ca9683f284d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4011704d-0c9f-4f7b-b1c8-ccd4898b2c3d",
   "metadata": {},
   "source": [
    "The first row in the results table is a parent source with 13 deblended children (it has `detect_isPrimary` = False, and `deblend_nChild` = 13).\n",
    "\n",
    "Double check that all of the children are included in the `table` by finding out how many rows have a `parentObjectId` matching the first row's `objectId`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0353e6b-7217-488f-b0a4-1f5ad408d4b9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sel_objid = table['objectId'][0]\n",
    "whchild = np.where(table['parentObjectId'] == sel_objid)[0]\n",
    "print(len(whchild), \"child objects deblended from parentObjectId\", sel_objid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "555bd9c2-2281-4eaf-b132-5809271b55b6",
   "metadata": {},
   "source": [
    "Up until now, the TAP service has been used to access catalog data.\n",
    "\n",
    "Below, all objects in the same tract and patch as the 13-child parent are obtained via the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f06ce5c0-d56f-4879-a5eb-b4cb5508d3cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tract = table['tract'][0]\n",
    "patch = table['patch'][0]\n",
    "dataId = {'tract': tract, 'patch': patch, 'band': 'i'}\n",
    "obj = butler.get('objectTable', dataId)\n",
    "print(\"Retrieved catalog of \", len(obj), \" objects.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8feb25e8-1730-4809-95e9-596fd7e1576b",
   "metadata": {},
   "source": [
    "Print the parameters only for the objects that are children of the parent identified above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74ae2755-16ea-4457-99b5-e12f8f6091ad",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "children = obj[obj.parentObjectId == sel_objid]\n",
    "children[['parentObjectId', 'x', 'y', 'coord_ra', 'coord_dec']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db481884-aa79-4bf5-b1c6-cdba308e3670",
   "metadata": {},
   "source": [
    "An alternative way of printing the `children` array defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03313095-8610-46f4-9efb-ef8cdccd8a82",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# for i in range(len(children)):\n",
    "#     print(children.iloc[i].parentObjectId, children.index[i],\n",
    "#           children.iloc[i].x, children.iloc[i].y,\n",
    "#           children.iloc[i].coord_ra, children.iloc[i].coord_dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd72077-ab48-4cc1-97db-69076f9a9cb7",
   "metadata": {},
   "source": [
    "Another alternative for defining separate arrays containing the data from `children`,\n",
    "and printing them with a user-specified format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2594ddec-33e9-41b1-9df5-f0f74277c4b0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# obj_parentObjectId = np.asarray(obj.parentObjectId, dtype='int')\n",
    "# obj_objectId = np.asarray(obj.index, dtype='int')\n",
    "# obj_x = np.asarray(obj.x, dtype='float')\n",
    "# obj_y = np.asarray(obj.y, dtype='float')\n",
    "# obj_ra = np.asarray(obj.coord_ra, dtype='float')\n",
    "# obj_dec = np.asarray(obj.coord_dec, dtype='float')\n",
    "\n",
    "# wh = np.where(obj_parentObjectId == sel_objid)[0]\n",
    "# print(len(wh), \"child objects found for parentObjectId\", sel_objid)\n",
    "\n",
    "# for i in wh:\n",
    "#     print('%20s %20s %10.3f %9.3f %10.6f %11.6f' %\n",
    "#           (obj_parentObjectId[i], obj_objectId[i], obj_x[i], obj_y[i],\n",
    "#            obj_ra[i], obj_dec[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35f52616-8403-47b2-84ff-729ba755f8f0",
   "metadata": {},
   "source": [
    "Obtain a cutout at the user-specified location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "572ed22d-d0bc-4f91-bae4-4778ac191d4f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cutout = cutout_coadd(butler, ra, dec, band='i', datasetType='deepCoadd',\n",
    "                      cutoutSideLength=201)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8270c096-ca7c-43de-899e-b062065ce822",
   "metadata": {},
   "source": [
    "Display a cutout at the location of this parent, and use the object data retrieved from the butler to mark the child locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31118cde-b512-4286-bade-4e1cd3e2df16",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout.image)\n",
    "with display.Buffering():\n",
    "    display.dot('o', table['x'][0], table['y'][0], size=5,\n",
    "                ctype=afwDisplay.YELLOW)\n",
    "    for i in range(len(children)):\n",
    "        display.dot('o', children.iloc[i].x, children.iloc[i].y,\n",
    "                    size=6, ctype=afwDisplay.GREEN)\n",
    "plt.gca().axis('off')\n",
    "display.show_colorbar(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec5d3f8-a5b4-4b4a-84e6-5f71a6fa833e",
   "metadata": {},
   "source": [
    "> Figure 2: A cutout image displayed in grayscale, with a yellow circle to mark the parent location and green circles to mark the children."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ea51c11-ad87-40b4-99b5-4c29130e2de1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del table, obj, cutout"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "086660d6-5013-4186-acaa-d7bd333be923",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. Deblended footprints\n",
    "\n",
    "Scarlet, when run as part of the LSST Science Pipelines, produces the `deepCoadd_deblendedFlux` data product (among other outputs) after using the multi-wavelength images to identify peaks.\n",
    "This data product contains the footprints for the individual deblended children, which might be a useful data product for some users.\n",
    "\n",
    "The concept of a footprint was introduced in tutorial 05.\n",
    "To quote [Bosch et al. (2018)](https://ui.adsabs.harvard.edu/abs/2018PASJ...70S...5B/abstract), \n",
    "\n",
    "> *Footprints record the exact above-threshold detection region on a CCD. These are similar to  Source Extractor’s “segmentation map\", in that they identify which pixels belong to which detected objects.*\n",
    "\n",
    "This quote draws an analogy between footprints and segmentation maps, since they both identify pixels with values above some threshold. This is a useful similarity, since it gives us a place to start understanding the properties of footprints.\n",
    "\n",
    "\n",
    "### 3.1. The `deepCoadd_deblendedFlux` table\n",
    "\n",
    "Below, find the `deepCoadd_deblendedFlux` data product in the butler by searching for all dataset types with the word \"blend\" in the title."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81d1d578-7c2a-433d-9d89-3d56f76c2940",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "registry = butler.registry\n",
    "for dt in sorted(registry.queryDatasetTypes('*blend*')):\n",
    "    print(dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3132683-9ed4-459a-9e50-19641e94d5d3",
   "metadata": {},
   "source": [
    "Look at the schema for the `deepCoadd_deblendedFlux` catalog to learn more about its contents.\n",
    "\n",
    "Use the `dataId` as defined above for the 13-child parent galaxy of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f74a74a9-d9af-4228-bcb4-1e7a82a218c5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dataId"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2e161b3-536c-4b83-bcac-5e3b9eb12e2d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-03-07T00:08:07.430297Z",
     "iopub.status.busy": "2023-03-07T00:08:07.430113Z",
     "iopub.status.idle": "2023-03-07T00:08:07.959780Z",
     "shell.execute_reply": "2023-03-07T00:08:07.957255Z",
     "shell.execute_reply.started": "2023-03-07T00:08:07.430282Z"
    }
   },
   "source": [
    "We can first look at the contexts of `deepCoadd_deblendedFlux` by looking at the column names in `deepCoadd_deblendedFlux_schema` without having to pull the contents for the dataID from the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e8af3d5-a88f-49cc-97f6-5b52feba2769",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dbFlux_schema = butler.get('deepCoadd_deblendedFlux_schema', dataId=dataId)\n",
    "for colname in dbFlux_schema.asAstropy().colnames:\n",
    "    print(colname)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb87dd4-063a-47dd-9ab6-163e63a01327",
   "metadata": {},
   "source": [
    "Obtain the catalog contents for the `deepCoadd_deblendedFlux` table in the tract and patch defined by our `dataId`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6657c83-ddb5-4a47-99ed-c393c686d58f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dbFlux = butler.get('deepCoadd_deblendedFlux', dataId=dataId)\n",
    "dbFlux_table = dbFlux.asAstropy()\n",
    "print('Retrieved', len(dbFlux_table),\n",
    "      'records from the deepCoadd_deblendedFlux table.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2ceee02-3d56-44a5-8e29-01bfdfa6fdc6",
   "metadata": {},
   "source": [
    "For the 13 children associated with the parent identified in Section 2.2, display the catalog content."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcb3660c-a14f-4cba-b19b-fc69f9dad62e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = np.where(dbFlux_table['parent'] == sel_objid)[0]\n",
    "dbFlux_table[tx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf89ba61-6e4f-4087-ad94-543c742903ba",
   "metadata": {},
   "source": [
    "### 3.2. Exploring deblended footprints\n",
    "\n",
    "Use the `.getFootprint()` method to get the footprints for all records in `dbFlux`,\n",
    "and add them to the astropy table made from the `deepCoadd_deblendedFlux` table.\n",
    "A new column called 'footprint' will appear in the `dbFlux_table`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2f5d229-71d9-4040-80ad-16f6d4c38dfb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "footprints = [src.getFootprint() for src in dbFlux]\n",
    "dbFlux_table[\"footprint\"] = footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f2c80bb-d6ad-4fa4-bb21-066b69b4ddf1",
   "metadata": {},
   "source": [
    "For the 13 children associated with the parent identified in Section 2.2, display the catalog content, which now contains a `footprint` column (scroll all the way over to the right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7523c2cc-6460-45fd-8a2c-c4c02452e8c1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = np.where(dbFlux_table['parent'] == sel_objid)[0]\n",
    "dbFlux_table[tx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3badb663-5341-4582-a2c2-0de8108c3b83",
   "metadata": {},
   "source": [
    "Use the function `heavyFootprint2Image` to turn the `footprint` into a `heavyFootprint` image. \n",
    "The function returns an image defined by a bounding box, and padded with fill value (the default is `NaN`)\n",
    "which contains the projected flux in `heavy`.\n",
    "\n",
    "The term \"heavy\" here refers to footprint pixels that aren't simply populated with values 0 or 1\n",
    "(where 1 indicates the pixel is part of the footprint),\n",
    "but with values that represent the flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83e06c7f-5827-41d0-b006-051c6b8c084f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fpimg = heavyFootprint2Image(dbFlux_table['footprint'][tx[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b54b6bb4-6305-40ab-a211-d56161974a77",
   "metadata": {},
   "source": [
    "*Optional:* Explore the properties of `fpimg` by uncommenting the line below, putting your cursor after the period, pressing tab, and scrolling through the pop-up box of attributes.\n",
    "Recognize familiar image attributes such as `image` and `getBBox`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbd476e2-d4fc-40ba-b653-2b5666c6d4eb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# fpimg."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22a784f4-411f-4f28-a4cb-b671f5774d08",
   "metadata": {},
   "source": [
    "Display the footprint of each of the 13 deblended children in a 3x5 grid.\n",
    "\n",
    "Title the subplot with the index and the footprint shape.\n",
    "Notice that the footprint image areas are not a standard size: they encompass the footprint.\n",
    "\n",
    "Remove the axis labels (x and y pixels) and the colorbar scale.\n",
    "Note that this grid leaves 2 subplots empty, so erase the axes for the final two subplots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6a53695-0f9e-43c9-bb3a-11bb609c526b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(3, 5, figsize=(10, 7))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "\n",
    "x = 0\n",
    "for i in range(3):\n",
    "    if i < 2:\n",
    "        Nc = 5\n",
    "    elif i == 2:\n",
    "        Nc = 3\n",
    "    for j in range(Nc):\n",
    "        fpimg = heavyFootprint2Image(dbFlux_table['footprint'][tx[x]])\n",
    "        plt.sca(ax[i, j])\n",
    "        display.scale('asinh', 'zscale')\n",
    "        display.mtv(fpimg.image)\n",
    "        display.show_colorbar(False)\n",
    "        plt.title(str(x) + \" size:\" + str(np.shape(fpimg.image.array)),\n",
    "                  fontsize=10)\n",
    "        plt.gca().axis('off')\n",
    "        x += 1\n",
    "        del fpimg\n",
    "\n",
    "plt.sca(ax[2, 3])\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "plt.gca().axis('off')\n",
    "plt.sca(ax[2, 4])\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "plt.gca().axis('off')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000a58ae-333a-429f-8ec8-095905075f8b",
   "metadata": {},
   "source": [
    "> Figure 3: The heavy footprints, displayed in grayscale, for each of the 13 children.\n",
    "\n",
    "The above grid is somewhat useful for seeing the shapes of the individual footprints, but it does not show how the footprints might overlap with each other.\n",
    "\n",
    "Below, plot all of the footprints together in a single axis. \n",
    "\n",
    "Define a set of 13 colormaps to use, one for each child footprint, and label each child footprint 0 through 12.\n",
    "\n",
    "The `alpha` parameter is used here to set the pixel transparency in a fairly crude way, so that pixels are semi-transparent in the bright (high flux) areas of the footprint, and totally transparent in the outskirts of the footprint.\n",
    "\n",
    "To make sure the final plot is big enough to hold all footprints, keep track of the min/max x/y values, and then use them to set the plot limits.\n",
    "\n",
    "> **Note:** Due to the scaling, it *seems* like not all of these footprints touch, but since they are all deblended from a single parent, they do in fact join and overlap.\n",
    "Adjust the values of `alpha` or try different `vmin` and `vmax` scaling to explore the footprint edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3303043-d458-4cc2-99a4-2f13dd0eb0f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_colormaps = ['Greys', 'Purples', 'Blues', 'Greens', 'Oranges',\n",
    "                'Reds', 'YlOrBr', 'RdPu', 'GnBu', 'YlGnBu',\n",
    "                'viridis', 'cividis', 'cool']\n",
    "xmin = 9999999.\n",
    "xmax = 0.\n",
    "ymin = 9999999.\n",
    "ymax = 0.\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "for x in range(len(tx)):\n",
    "    fpimg = heavyFootprint2Image(dbFlux_table['footprint'][tx[x]])\n",
    "    fpimg_extent = (fpimg.getBBox().beginX, fpimg.getBBox().endX,\n",
    "                    fpimg.getBBox().beginY, fpimg.getBBox().endY)\n",
    "    alphas = np.asarray(fpimg.image.array,\n",
    "                        dtype='float') / np.nanmax(fpimg.image.array)\n",
    "    tmp_x = np.where((np.isnan(alphas)) | (alphas < 0.1))\n",
    "    alphas[tmp_x] = 0.0\n",
    "    tmp_x = np.where(alphas >= 0.1)\n",
    "    alphas[tmp_x] = 0.7\n",
    "    im = plt.imshow(fpimg.image.array, cmap=my_colormaps[x], alpha=alphas,\n",
    "                    vmin=np.nanmin(fpimg.image.array),\n",
    "                    vmax=np.nanmax(fpimg.image.array),\n",
    "                    extent=fpimg_extent, origin='lower')\n",
    "    if xmin > fpimg.getBBox().beginX:\n",
    "        xmin = fpimg.getBBox().beginX\n",
    "    if xmax < fpimg.getBBox().endX:\n",
    "        xmax = fpimg.getBBox().endX\n",
    "    if ymin > fpimg.getBBox().beginY:\n",
    "        ymin = fpimg.getBBox().beginY\n",
    "    if ymax < fpimg.getBBox().endY:\n",
    "        ymax = fpimg.getBBox().endY\n",
    "    del fpimg, fpimg_extent, alphas, tmp_x, im\n",
    "\n",
    "    center_x = dbFlux_table['deblend_peak_center_x'][tx[x]]\n",
    "    center_y = dbFlux_table['deblend_peak_center_y'][tx[x]]\n",
    "    plt.text(center_x, center_y, str(x), fontsize=10)\n",
    "\n",
    "plt.xlim([xmin, xmax])\n",
    "plt.ylim([ymin, ymax])\n",
    "plt.xlabel('x', fontsize=10)\n",
    "plt.ylabel('y', fontsize=10)\n",
    "plt.tick_params(labelsize=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16424d9a-5841-408d-b604-1001a9ca0861",
   "metadata": {},
   "source": [
    "> Figure 4: The heavy footprints from Figure 3, each with their own colormap applied (instead of grayscale) and labeled as 0 through 12, to distinguish them. Instead of showing each heavy footprint separately, they've been pasted back into the image frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c47315-cdd2-4ffa-9ee6-642981199166",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del dataId, dbFlux_schema, dbFlux, dbFlux_table\n",
    "del footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea7684a5-df94-46b1-8780-9ac09a02046b",
   "metadata": {},
   "source": [
    "## 4. Use-case: SN Ia host galaxy identification\n",
    "\n",
    "DP0 <a href=\"https://dp0-2.lsst.io/tutorials-examples/portal-images.html\">Portal tutorial 03</a> \n",
    "investigated the location of a Type Ia supernova (SN Ia) near a bright host galaxy and found that in some filters,\n",
    "and with some image scaling parameters, it seemed like the true host is a fainter host\n",
    "that is blended with a bigger, brighter galaxy.\n",
    "\n",
    "The following steps use the techniques above to look for parent and child sources near this SN Ia.\n",
    "\n",
    "Define the coordinates for this SN Ia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d99ea6-5a45-4c6e-8954-5f9bdb5797d5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sn_ra = 67.4579\n",
    "sn_dec = -44.0802"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "501409bf-ba44-4cb9-9a08-c401e3376aad",
   "metadata": {},
   "source": [
    "Obtain the i-band `deepCoadd` image, which in Step 2.3 of Portal tutorial 03 reveals a potentially blended host."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9897101-70b6-48f9-9690-b5351c6a3911",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cutout = cutout_coadd(butler, sn_ra, sn_dec, band='i', datasetType='deepCoadd',\n",
    "                      cutoutSideLength=71)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe7c38c1-0c35-4700-8140-24263a0bb955",
   "metadata": {},
   "source": [
    "Get the WCS for this cutout and convert the SN Ia coordinates from RA, Dec to pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b5918e2-7fda-4932-b90f-6024662795c1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wcs = cutout.getWcs()\n",
    "sn_radec = geom.SpherePoint(sn_ra, sn_dec, geom.degrees)\n",
    "sn_xy = wcs.skyToPixel(sn_radec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ca9e404-d453-4dd9-a1ef-3d4af0ad0a7c",
   "metadata": {},
   "source": [
    "Display the i-band `deepCoadd` with the SN Ia location marked as a pink circle.\n",
    "\n",
    "See how it seems like there's a faint galaxy there, blended with the larger one?\n",
    "\n",
    "Try changing the scale from `asinh` to `linear` to get a different view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fec17d7-d65c-4e13-9791-182b9e70e5c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 4))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout.image)\n",
    "with display.Buffering():\n",
    "    display.dot('o', sn_xy.getX(), sn_xy.getY(), size=4,\n",
    "                ctype=afwDisplay.MAGENTA)\n",
    "\n",
    "plt.gca().axis('off')\n",
    "display.show_colorbar(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7414d8b0-a0ab-4766-b33e-7c78891ba0e9",
   "metadata": {},
   "source": [
    "> Figure 5: The cutout, in grayscale, at the location of a supernova (magenta circle).\n",
    "\n",
    "Query the `Object` catalog for extended parents in the area. As a reminder, the spatial constraint by searching around a specified ra/dec is not technically necessary, but speeds up the query significantly because of the way the data in the object Table are stored (see Section 2.1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3369cf9f-07aa-4c65-966d-50673527e8ae",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, coord_ra, coord_dec, x, y, tract, \" + \\\n",
    "        \"patch, refExtendedness, detect_isPrimary, detect_isIsolated, \" + \\\n",
    "        \"deblend_nChild, detect_fromBlend, footprintArea, parentObjectId \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec ), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + str(sn_ra) + \", \" + str(sn_dec) + \\\n",
    "        \", 0.005)) = 1 \" + \\\n",
    "        \"AND refExtendedness = 1 AND deblend_nChild > 0\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eca88fef-a66a-4e49-be5f-17caa39f2f5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02cb114a-ad5c-449d-975b-40f3ab693da8",
   "metadata": {},
   "outputs": [],
   "source": [
    "parents = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8286ea9-6aea-482e-b964-64e005aed1f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "parents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16ae2c4b-7614-4270-aaa6-8fa2df0c4e34",
   "metadata": {},
   "source": [
    "The table above reveals there is one nearby parent, with 32 children (identified by the deblend_nChild keyword). The other 5 sources you can see are not blends because they have detect_isIsolated = True.\n",
    "\n",
    "Plot the location of this parent on the `deepCoadd` to see if it's the bright nearby galaxy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7730bcd-5a0e-49e1-a110-a9983c316442",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = np.where(parents['deblend_nChild'] == 32)[0]\n",
    "print('objectId of parent with 32 children: ', parents['objectId'][tx[0]])\n",
    "\n",
    "fig = plt.figure(figsize=(4, 4))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout.image)\n",
    "with display.Buffering():\n",
    "    display.dot('o', sn_xy.getX(), sn_xy.getY(), size=4,\n",
    "                ctype=afwDisplay.MAGENTA)\n",
    "    display.dot('o', parents['x'][tx], parents['y'][tx], size=5,\n",
    "                ctype=afwDisplay.YELLOW)\n",
    "\n",
    "plt.gca().axis('off')\n",
    "display.show_colorbar(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8be4f10b-eb88-4de2-9e27-3f4e67e5a806",
   "metadata": {},
   "source": [
    "> Figure 6: Same as Figure 5, but with a yellow circle to mark a nearby parent object.\n",
    "\n",
    "Query for the 32 children of this parent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b122397b-c985-4132-b431-d7472ba27237",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, coord_ra, coord_dec, x, y, tract, patch, \" + \\\n",
    "        \"refExtendedness, detect_isPrimary, deblend_nChild, \" + \\\n",
    "        \"detect_fromBlend, footprintArea, parentObjectId \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec ), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + str(sn_ra) + \", \" + str(sn_dec) + \\\n",
    "        \", 0.01)) = 1 AND parentObjectId = 1252220598734559348\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7a2fac1-0f22-4597-b03b-e49217bf2a7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "200b6d96-4c96-4690-a45a-03b8420b2161",
   "metadata": {},
   "outputs": [],
   "source": [
    "children = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44401f74-5429-43ea-a6a2-093b332b90e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# children"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0361426-bc27-4220-bddd-9767403df76c",
   "metadata": {},
   "source": [
    "Plot the locations of the children as green circles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0cb5f4c-56ad-48bd-9038-c944e0c0b0b4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 4))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(cutout.image)\n",
    "with display.Buffering():\n",
    "    display.dot('o', sn_xy.getX(), sn_xy.getY(), size=4,\n",
    "                ctype=afwDisplay.MAGENTA)\n",
    "    display.dot('o', parents['x'][tx], parents['y'][tx], size=5,\n",
    "                ctype=afwDisplay.YELLOW)\n",
    "    for i in range(len(children)):\n",
    "        display.dot('o', children['x'][i], children['y'][i], size=3,\n",
    "                    ctype=afwDisplay.GREEN)\n",
    "\n",
    "plt.gca().axis('off')\n",
    "display.show_colorbar(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5527e6c6-d1a1-4f31-8a80-d6f969cf57df",
   "metadata": {},
   "source": [
    "> Figure 7: Same as Figure 6, but with green circles to mark the children.\n",
    "\n",
    "Well! What seems to be the true host galaxy has not been deblended from the large, bright, nearby galaxy.\n",
    "\n",
    "This might be a good case where the deblender should be run, with more aggressive deblending parameters,\n",
    "but that is an intermediate/advanced level use of the LSST Science Pipelines that is left to a different tutorial.\n",
    "\n",
    "## 5. Exercise for the learner\n",
    "\n",
    "Further explore the footprints of the children around this SN Ia.\n",
    "Create a footprints plot similar to the one in Section 3.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1368ea07-b7d7-4447-a39a-aeedb5c61509",
   "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.9"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}