DP02_05_Source_Detection_and_Measurement.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5462d256-8c6c-419e-a5a9-3213e9d28ced",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-17T17:16:09.765000Z",
     "iopub.status.busy": "2024-05-17T17:16:09.764040Z",
     "iopub.status.idle": "2024-05-17T17:16:09.770726Z",
     "shell.execute_reply": "2024-05-17T17:16:09.769774Z",
     "shell.execute_reply.started": "2024-05-17T17:16:09.764966Z"
    }
   },
   "source": [
    "# Source Detection and Measurement\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>Contact authors: Douglas Tucker, Alex Drlica-Wagner<br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_50 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: Intermediate <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** Access, display, and manipulate images; detect, deblend, and measure sources; and extract, plot, and use object footprints."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Run source detection, deblending, and measurement. Use source footprints."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** DP0.2 processed visit images and catalogs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** lsst.pipe.tasks.characterizeImage, lsst.meas.algorithms.detection.SourceDetectionTask, lsst.meas.deblender.SourceDeblendTask, lsst.meas.base.SingleFrameMeasurementTask, lsst.meas.base.ForcedMeasurementTask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:**\n",
    "Originally developed by Alex Drlica-Wagner and Imran Hasan in the context of the LSST Stack Club.\n",
    "Forced photometry routine by Tobias Géron and Melissa Graham "
   ]
  },
  {
   "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 provides a brief introduction to running the LSST Science Pipelines source detection, deblending, and measurement tasks.\n",
    "It does not go into depth about optimizing these tasks for different types of sources.\n",
    "It also covers how to access the footprint of a detected source and how to do forced photometry at user-specified coordinates.\n",
    "\n",
    "Some source detection and measurement details come from Robert Lupton's [Tune Detection.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Tune%20Detection.ipynb) and [Kron.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Kron.ipynb).\n",
    "Interaction with `lsst.afw.display` was also improved by studying Michael Wood-Vasey's [DC2_Postage Stamps.ipynb](https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/dm_butler_postage_stamps.ipynb).\n",
    "More information on footprints can be found in [this Stack Club notebook by Imran Hasan](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/SourceDetection/Footprints.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1. Import packages\n",
    "\n",
    "The [`numpy`](https://numpy.org/) and [`matplotlib`](https://matplotlib.org/) are widely used Python plot libraries.\n",
    "\n",
    "From the [LSST Science Pipelines](https://pipelines.lsst.io/), modules for data access via the butler, image display,\n",
    "image geometry, table handling, and various tasks for image characterization and source detection, deblending,\n",
    "and measurement are imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc506441-c5c1-4101-9090-3f2c59a04a37",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import getpass\n",
    "\n",
    "import lsst.daf.base as dafBase\n",
    "from lsst.daf.butler import Butler\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.afw.table as afwTable\n",
    "import lsst.geom as geom\n",
    "\n",
    "from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask\n",
    "from lsst.meas.algorithms.detection import SourceDetectionTask\n",
    "from lsst.meas.deblender import SourceDeblendTask\n",
    "from lsst.meas.base import SingleFrameMeasurementTask\n",
    "from lsst.meas.base import ForcedMeasurementTask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce4e490-ad26-4738-b555-920096deec40",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e6e035b-e492-42ed-b57a-ddbcff9a3220",
   "metadata": {},
   "source": [
    "Set `afwDisplay` to use `Firefly` as its backend.\n",
    "\n",
    "Set `afw_display` to be frame 1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "747c7f53-2418-4784-ae3b-375e61b07a52",
   "metadata": {},
   "source": [
    "Firefly should have opened in a new tab in the same workspace area as this notebook.\n",
    "\n",
    "Click on the \"Firefly: stateClient...\" tab and drag it down and to the right of this workspace in order to work with the notebook and Firefly side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36096062-7851-4099-bd5a-ad35efacb434",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('firefly')\n",
    "afw_display = afwDisplay.Display(frame=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bd0bd8f-d8aa-40e4-bd17-8c38b8b00afe",
   "metadata": {},
   "source": [
    "Set the `matplotlib` plot color table to be colorblind-friendly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d36ac471-cf85-4621-b89d-2107617d66c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd072d0c-2136-426c-b934-62f914266fe1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-01T00:49:45.669777Z",
     "iopub.status.busy": "2024-06-01T00:49:45.669628Z",
     "iopub.status.idle": "2024-06-01T00:49:45.700618Z",
     "shell.execute_reply": "2024-06-01T00:49:45.700211Z",
     "shell.execute_reply.started": "2024-06-01T00:49:45.669764Z"
    },
    "tags": []
   },
   "source": [
    "Instantiate butler access to the DP0.2 data sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad1312ca-a401-4fa2-ba97-c4d271ab8c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "butler = Butler('dp02', collections='2.2i/runs/DP0.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a26304af-d7a5-4625-8adf-294612a7b654",
   "metadata": {},
   "source": [
    "## 2. Image access, display, and manipulation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5054c016-b2dc-491d-bac4-0a0ef1c6c87e",
   "metadata": {},
   "source": [
    "### 2.1. Retrieve and display a `calexp`\n",
    "\n",
    "Define the `dataId` for pre-selected visit 512055 and detector 75,\n",
    "then use the `butler` to retrieve the `calexp` (processed visit image)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3ee02e9-9951-44af-a3d9-8a32bf61b8fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dataId = {'visit': 512055, 'detector': 75}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6822c187-d6e3-43c6-a37e-90fe5bbcc2b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "calexp = butler.get('calexp', dataId=dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a164425-c575-4b86-9fa1-ad6095c67dc3",
   "metadata": {},
   "source": [
    "A `calexp` contains the image plus a mask and variance planes.\n",
    "\n",
    "**Option:** uncomment any of the following lines and execute the following cell to learn more about the various aspects of the `calexp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "599907bd-5328-40e2-8abc-57a521db8d9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calexp.maskedImage\n",
    "# calexp.maskedImage.image\n",
    "# calexp.maskedImage.mask\n",
    "# calexp.maskedImage.variance\n",
    "# calexp.image\n",
    "# calexp.mask\n",
    "# calexp.variance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edfc3e1c-82ea-44b3-8082-0cd9b39f06f1",
   "metadata": {},
   "source": [
    "The retrieved `calexp` also contains the point-spread function (PSF), the world coordinate system (WCS),\n",
    "and the photometric calibration for the image.\n",
    "\n",
    "These components could be retreived with the following statements.\n",
    "```\n",
    "psf = calexp.getPsf()\n",
    "wcs = calexp.getWcs()\n",
    "pcal = calexp.getPhotoCalib()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d634dfc5-90f7-4f17-b0c0-f5ffd52c9419",
   "metadata": {},
   "source": [
    "Display the `calexp` image and set the mask layer to be 100% transparent,\n",
    "so that the display shows only the image (i.e., does not color the pixels by mask value)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a988a6e0-735a-4b11-aa67-df4d5ca2577c",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display.mtv(calexp)\n",
    "afw_display.setMaskTransparency(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d581c53a-714f-4b54-a4c1-63bfc9efebb9",
   "metadata": {},
   "source": [
    "### 2.2. Add the sky background back into the image\n",
    "\n",
    "As an example of how to manipulate an image before reprocessing, use the `butler` to\n",
    "retrieve the sky background which was subtracted, and add it back in to the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c96653ac-20ec-4343-a195-63b331cc6713",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bkgd = butler.get('calexpBackground', dataId=dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cb47b27-e037-49e9-968e-6f5d52e3c6e4",
   "metadata": {},
   "source": [
    "Reset `afw_display` to be for frame 2, and display the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f18435-6037-4cbb-90d8-36460b1ab11a",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display = afwDisplay.Display(frame=2)\n",
    "afw_display.mtv(bkgd.getImage())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebca7a86-c50d-407a-bf12-f6e18bfb0be4",
   "metadata": {},
   "source": [
    "Add the background into the `calexp` and display the results.\n",
    "\n",
    "**Warning!** Only execute this cell once or the background will be re-added multiple times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bdd474d-d6c8-48ad-99f6-f662b8cc94cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "calexp.maskedImage += bkgd.getImage()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0a31e1f-73b5-4abe-9464-6fb82f06de1c",
   "metadata": {},
   "source": [
    "Reset `afw_display` to be for frame 3, display the `calexp`,\n",
    "and set mask transparency to 100%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47d1a3b9-33f6-44f4-840c-d77ead5b0d7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display = afwDisplay.Display(frame=3)\n",
    "afw_display.mtv(calexp)\n",
    "afw_display.setMaskTransparency(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6937678e-2536-49b0-8773-067995f10adc",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/dp02_05_fig1.png\" width=500>\n",
    "\n",
    "> Figure 1: A screenshot of what the Firefly interface should look like with the original `calexp` (upper left), the background image (upper right), and the `calexp` with the background added back in (lower left) all displayed in frames 1 through 3, respectively.\n",
    "\n",
    "**If images are displayed side-by-side instead of tiled**,\n",
    "as in the screenshot above,\n",
    "find the \"Scroll Images\" toggle at upper left in the Firefly window and click it.\n",
    "\n",
    "Compare the pixel values of the original `calexp` in frame 1 and the\n",
    "new `calexp` with background in frame 3.\n",
    "\n",
    "**Action:** \n",
    "With the mouse, hover over any image in the Firefly window and see the coordinates\n",
    "and pixel value for the mouse position appear at the bottom of the image.\n",
    "\n",
    "If the pixel value does not appear, only coordinates, there are two main options.\n",
    "\n",
    "1. Expand the size of the Firefly window by clicking and dragging the grey vertical line between tutorial notebook and Firefly window until \"Value\" appears to the right of \"EQ-J2000\".\n",
    "\n",
    "2. Click on the \"Show expanded readout\" icon in the lower-left corner (a box with an arrow pointed at the upper-right corner). A pop-up window will open which contains the pixel value at the bottom.\n",
    "\n",
    "Now with the mouse, hover over the images in frames 1, 2, and 3 and see how the pixel values differ by the value of the background pixels.\n",
    "\n",
    "Reset the Firefly window size or close the pop-up window.\n",
    "\n",
    "**Action:** Close each frame by clicking on the X in the upper-right corner of the image (it is not very visible until the mouse is over it). The Firefly interface should now say \"Firefly Ready\"."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "603fb547-b971-402b-8cf1-7d73c852afa2",
   "metadata": {},
   "source": [
    "Clean up. Do not use the `calexp` + `bkgd` in the sections below. A fresh version of the `calexp` is retrieved in Section 3.1.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0490c4a-db2e-4efb-a81c-b4adec0fda75",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del calexp, bkgd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a7d475f-df1f-476b-89cf-9f912833c31a",
   "metadata": {},
   "source": [
    "## 3. Source detection, deblending, and measurement"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caa35bd1-6e37-4aa9-b78d-87f958dee916",
   "metadata": {},
   "source": [
    "The LSST Science Pipelines' source detection, deblending, and measurement tasks were imported in the first code cell.\n",
    "\n",
    "More information can be found at [pipelines.lsst.io](https://pipelines.lsst.io/), either by scrolling down and\n",
    "browsing the available modules or by using the search bar at upper left.\n",
    "\n",
    "The schema describes the output properties that will be measured for each source.\n",
    "The schema needs to be passed to all of the tasks, each of which will add columns to it when it runs.\n",
    "\n",
    "Create a minimal schema for a source table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aa6787b-aaa5-44a3-ba6c-db3ad87a7bf4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "schema = afwTable.SourceTable.makeMinimalSchema()\n",
    "print(schema)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4b93ace-5162-4dda-a1b1-dbea9ab80d6a",
   "metadata": {},
   "source": [
    "Add `coord_raErr` and `coord_decErr` to the minimal schema to avoid warnings \n",
    "from `sourceMeasurementTask` in Section 3.4, which expects these columns to exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7b44bfa-948b-467f-8555-0079e2a0bc32",
   "metadata": {},
   "outputs": [],
   "source": [
    "raerr = schema.addField(\"coord_raErr\", type=\"F\")\n",
    "decerr = schema.addField(\"coord_decErr\", type=\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3adcc13-8ae8-44b1-ad5d-6c3b1bf0f004",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-02T00:17:55.647565Z",
     "iopub.status.busy": "2024-06-02T00:17:55.647343Z",
     "iopub.status.idle": "2024-06-02T00:17:55.650905Z",
     "shell.execute_reply": "2024-06-02T00:17:55.650325Z",
     "shell.execute_reply.started": "2024-06-02T00:17:55.647550Z"
    }
   },
   "source": [
    "Create a container which will be used to record metadata about algorithm execution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d6b1402-bd72-411b-975d-423a2b2f982a",
   "metadata": {},
   "outputs": [],
   "source": [
    "algMetadata = dafBase.PropertyList()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cf0992b-f964-40ff-8f65-1f1b1a6de733",
   "metadata": {},
   "source": [
    "### 3.1. Configure the tasks\n",
    "\n",
    "Each task possesses an associated configuration class. The properties of these configuration classes can be determined from the classes themselves.\n",
    "\n",
    "**Option:** uncomment the following line to view the help documentation for the `CharacterizeImageTask` configuration.\n",
    "Replace `CharacterizeImageTask` with another task name to view its help documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8603996-260e-4dd9-88f5-b1c74d1ff3e9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# help(CharacterizeImageTask.ConfigClass())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76f84d34-43aa-42d6-a5dc-0e9e5359d5ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# help(SourceDetectionTask.ConfigClass())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddac4a8e-78d3-48c2-a6ea-3c05239c92b5",
   "metadata": {},
   "source": [
    "Set the basic configuration parameters and instantiate each task in turn.\n",
    "\n",
    "`CharacterizeImageTask` characterizes the image properties (e.g., the PSF).\n",
    "\n",
    "Set the configuration to only do one PSF iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bced8cd2-a326-4c96-aeae-bdc31bbc3bcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = CharacterizeImageTask.ConfigClass()\n",
    "config.psfIterations = 1\n",
    "charImageTask = CharacterizeImageTask(config=config)\n",
    "del config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "675736b4-aa6e-4f51-b775-2580b70732c6",
   "metadata": {},
   "source": [
    "`SourceDetectionTask` detects sources.\n",
    "\n",
    "Set the configuration for the detection threshold to 4 (the default for LSST catalogs is 5) in order to detect sub-threshold (faint) sources.\n",
    "Set the detection threshold type to be `stdev`, which means the threshold is applied to the image's standard deviation.\n",
    "In other words, this configuration will detect sources at $4\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb17104c-aeb0-46fd-894d-cbdccedf8b17",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = SourceDetectionTask.ConfigClass()\n",
    "config.thresholdValue = 4\n",
    "config.thresholdType = \"stdev\"\n",
    "sourceDetectionTask = SourceDetectionTask(schema=schema, config=config)\n",
    "del config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8d61ce9-c821-4edb-a26f-4336eb8647f9",
   "metadata": {},
   "source": [
    "`SourceDeblendTask` deblends overlapping sources into \"parent\" and \"child\" sources.\n",
    "No configuration is needed for this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93e95ad6-1a52-4d35-b112-3d226e95afdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "sourceDeblendTask = SourceDeblendTask(schema=schema)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f31995a6-7972-4de8-81f0-b15b10ba3c85",
   "metadata": {},
   "source": [
    "`SingleFrameMeasurementTask` measures the properties of the deblended sources. No special configuration is needed for this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c434d66c-4ddc-42af-8b69-ee25a67a04ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "config = SingleFrameMeasurementTask.ConfigClass()\n",
    "sourceMeasurementTask = SingleFrameMeasurementTask(schema=schema,\n",
    "                                                   config=config,\n",
    "                                                   algMetadata=algMetadata)\n",
    "del config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "455b0e6c-2c15-4334-a402-9cffcf02ecce",
   "metadata": {},
   "source": [
    "Create a `SourceTable` called `tab` using the minimal schema. It will hold the output of the analysis in the following section.\n",
    "\n",
    "Then, as `schema` is no longer needed and a fresh minimal schema is made in Section 5, delete this one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2259f2ee-aa9c-4109-bf66-63d7a3395a77",
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = afwTable.SourceTable.make(schema)\n",
    "del schema"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5a3f968-af7c-4f26-b1db-3f011ae87e2d",
   "metadata": {},
   "source": [
    "#### 3.1.1. Re-configure tasks and explore methods\n",
    "\n",
    "The configuration parameters cannot be changed after the task is constructed.\n",
    "To change a configuration parameter, redefine it and then also redefine the task.\n",
    "\n",
    "For example, to change the number of PSF iterations done during image characterization:\n",
    "```\n",
    "config.psfIterations = 3\n",
    "charImageTask = CharacterizeImageTask(config=config)\n",
    "```\n",
    "\n",
    "**Option:** To explore a task, use the `help` function on the task or any of its methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "625187c6-090d-4509-b487-5306e44ef667",
   "metadata": {},
   "outputs": [],
   "source": [
    "# help(charImageTask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a9ee579-7b51-46bc-8a21-eb5814fd54bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# help(charImageTask.writeSchemas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93939732-adef-45f8-b567-1d513df9e4c2",
   "metadata": {},
   "source": [
    "**Option:** To see a pop-up list of all methods for a given task, uncomment this line, place the cursor after the period, and press tab.\n",
    "Then re-comment the line. Do not execute the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1d18b04-3659-40d7-9dbd-974eec79091b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# charImageTask."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d349fb9-c61b-465e-b381-bc2798c76acf",
   "metadata": {},
   "source": [
    "#### 3.1.2. Retrieve a fresh version of the image\n",
    "\n",
    "Retrieve a fresh, unaltered `calexp` image from the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e1fc576-09f6-4991-a92c-3a8b853af845",
   "metadata": {},
   "outputs": [],
   "source": [
    "calexp = butler.get('calexp', dataId=dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dec41cee-04d9-47c6-beb5-78d8c731cfc9",
   "metadata": {},
   "source": [
    "Clear the \"DETECTED\" information from the mask plane, as source detection\n",
    "is going to be re-done below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9165704-11fb-4435-85a7-9c9707048e8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "calexp.mask.removeAndClearMaskPlane('DETECTED')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c767c92-e510-445d-bf52-11341ba0cfee",
   "metadata": {},
   "source": [
    "### 3.2. Image characterization\n",
    "\n",
    "Characterize the `calexp` image. This calculates various global properties of the image, such as the PSF FWHM (full-width half-max)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf0d483f-5876-4918-9b7b-c4c756699531",
   "metadata": {},
   "outputs": [],
   "source": [
    "result = charImageTask.run(calexp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff0df23-322a-4ffd-b500-392d002ecfc5",
   "metadata": {},
   "source": [
    "The `result` is a structure, [`lsst.pipe.base.struct.Struct`](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1pipe_1_1base_1_1struct_1_1_struct.html).\n",
    "\n",
    "A `Struct` is just a generalized container for storing multiple output components and accessing them as attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "574ada2a-7464-4a2e-8a20-88741b089461",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a29d305-4efe-4675-84f8-46517ffac387",
   "metadata": {},
   "source": [
    "The content of this `Struct` can be investigated with the `getDict` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68f631a-292d-4e5b-83d9-d5f3c6e34311",
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, v in result.getDict().items():\n",
    "    print(k, type(v))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df93c7b5-9470-45d5-b97d-75faa8de5976",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-02T01:08:41.690328Z",
     "iopub.status.busy": "2024-06-02T01:08:41.690114Z",
     "iopub.status.idle": "2024-06-02T01:08:41.693512Z",
     "shell.execute_reply": "2024-06-02T01:08:41.693060Z",
     "shell.execute_reply.started": "2024-06-02T01:08:41.690314Z"
    }
   },
   "source": [
    "**Option:** the `Struct` can be displayed but this is not very useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eea5850f-8ca4-4560-a157-40d2bd87981c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2aa50d7-37dd-4548-be09-19885fb6caa3",
   "metadata": {},
   "source": [
    "As an example of how to interact with the results of image characterization, get the PSF at a specific\n",
    "pixel location and print the FWHM. The factor of 2.355 converts from the standard deviation to the FWHM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef226bce-36cd-43ce-bb02-9cc9dcf19b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = 1700, 2100\n",
    "point = geom.Point2D(x, y)\n",
    "psf = calexp.getPsf()\n",
    "sigma = psf.computeShape(point).getDeterminantRadius()\n",
    "pixelScale = calexp.getWcs().getPixelScale().asArcseconds()\n",
    "print('PSF FWHM = {:.2f} arcsec'.format(sigma*pixelScale*2.355))\n",
    "\n",
    "del x, y, point, psf, sigma, pixelScale"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "071cb4bc-a5b0-4908-8fed-2242c71c5a1f",
   "metadata": {},
   "source": [
    "### 3.3. Source detection\n",
    "\n",
    "Run source detection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6566e107-559d-4eae-830d-e001e4fe42c1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "result = sourceDetectionTask.run(tab, calexp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "386d20f9-8234-4e0a-8bdf-c465fccb0b07",
   "metadata": {},
   "source": [
    "Print the number of positive peaks detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6ab3e88-432a-4f0e-9168-0c72956e969d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "result.numPosPeaks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89e9f774-001b-4558-9579-16723ab00e0d",
   "metadata": {},
   "source": [
    "With `thresholdValue` set to 4, over 2000 sources are detected.\n",
    "\n",
    "Had the `thresholdValue` been set instead to 10, there would have been about 940 peaks detected.\n",
    "\n",
    "Create `sources` to hold the detected sources from the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "444d2882-7758-4389-9d63-cda648e1920d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sources = result.sources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52c4d28d-a76a-4438-a07d-d3d3614828ed",
   "metadata": {},
   "source": [
    "**Option:** Display `sources` to see that most of the columns are `NaN` at first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c278486f-c8f0-45f4-9c2f-1b7ce54975d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# sources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e071a046-bb3c-4bf3-a667-14f1bbbd3996",
   "metadata": {},
   "source": [
    "**Option:** Write the `sources` to a FITS table, or save the `calexp` as a FITS file, in the home directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c872a0e9-e380-4833-93b5-5b4a61a67040",
   "metadata": {},
   "outputs": [],
   "source": [
    "# home_directory = '/home/' + getpass.getuser() + '/'\n",
    "# sources.writeFits(home_directory + \"DP02_NB05_sources.fits\")\n",
    "# calexp.writeFits(home_directory + \"DP02_NB05_calexp.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc659715-c636-43d8-b727-54d09936276f",
   "metadata": {},
   "source": [
    "### 3.4. Source debelending and measurement\n",
    "\n",
    "Run the `SourceDeblendTask` and `SingleFrameMeasurementTask`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4587f17-fef0-409d-8e68-e18a079ee94e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sourceDeblendTask.run(calexp, sources)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64e7fe05-773d-4ec8-8975-47d1f096ed6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sourceMeasurementTask.run(measCat=sources, exposure=calexp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e588432-ff0f-4fc7-bdcc-4d6cc47e26cd",
   "metadata": {},
   "source": [
    "Ensure that the `sources` catalog is contiguous (sequential) in memory by making a copy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d32012d-8bbd-4962-85c7-402b948c495e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sources = sources.copy(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f72f6074-7b41-4f63-a1af-05edeae6fd06",
   "metadata": {},
   "source": [
    "**Option:** Display the `sources` catalog as an `astropy` table, which has a nice, human-readable output format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "176682dd-a215-4fe9-88ab-bac7936e0357",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sources.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5760a2d-2eaf-45ba-818f-c7d50fe7011d",
   "metadata": {},
   "source": [
    "### 3.5. Visualize sources on an image cutout\n",
    "\n",
    "Overplot the detected sources on a cutout image of the `calexp` using `afwDisplay`.\n",
    "\n",
    "There are several ways to create image cutouts. This tutorial demonstrates the use of the `Factory` method for a `calexp`.\n",
    "\n",
    "The arguments that must be passed to the `Factory` method are:\n",
    "```\n",
    "calexp : the ExposureF to make the cutout from\n",
    "bbox   : the desired bounding box of the cutout\n",
    "origin : the image pixel origin, local to the cutout array\n",
    "deep   : if True, copy the data rather than passing by reference\n",
    "```\n",
    "\n",
    "Define the `bbox` region for a cutout that is centered on pixel coordinate 1700, 2100 and is 400 pixels wide and 400 pixels high."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64617ba5-304b-4fa8-90f9-2588bbcf19dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = 1700, 2100\n",
    "w, h = 400, 400\n",
    "xmin, ymin = x - w//2, y - h//2\n",
    "print(xmin, ymin)\n",
    "bbox = geom.Box2I()\n",
    "bbox.include(geom.Point2I(xmin, ymin))\n",
    "bbox.include(geom.Point2I(xmin + w, ymin + h))\n",
    "del x, y, w, h, xmin, ymin"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d771afc9-2fad-40a4-8633-854b3b142b51",
   "metadata": {},
   "source": [
    "**Option:** Display the `bbox` information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfafea79-e2ba-4fda-b84d-7d5d72efcc3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "870a299d-5c9c-4707-aa3b-f79d96863c44",
   "metadata": {},
   "source": [
    "An alternative statement to create `bbox` is:\n",
    "```\n",
    "bbox = geom.Box2I(geom.Point2I(xmin, ymin), geom.Extent2I(w, h))\n",
    "```\n",
    "\n",
    "Generate the cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d901b23-faf3-4d35-aa51-6fff33d02780",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cutout = calexp.Factory(calexp, bbox, origin=afwImage.LOCAL, deep=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf2636ab-d6d9-4c0c-b4fb-c1949f28b8be",
   "metadata": {},
   "source": [
    "Reset `afw_display` to be for frame 1, display the cutout, and\n",
    "set mask transparency to 100%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fdc10db-ba4a-4a9a-8812-9d029b8ff585",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display = afwDisplay.Display(frame=1)\n",
    "afw_display.mtv(cutout)\n",
    "afw_display.setMaskTransparency(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec10acce-3a0d-4870-947e-715f24875bf0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-10T20:52:53.259791Z",
     "iopub.status.busy": "2024-06-10T20:52:53.259042Z",
     "iopub.status.idle": "2024-06-10T20:52:53.322926Z",
     "shell.execute_reply": "2024-06-10T20:52:53.322329Z",
     "shell.execute_reply.started": "2024-06-10T20:52:53.259767Z"
    }
   },
   "source": [
    "If frame 1 only fills half of the Firefly interface, click and drag the lower right corner to resize frame 1.\n",
    "\n",
    "Overplot sources as orange circles.\n",
    "Use display buffering while overplotting sources to avoid re-drawing the image after each source is plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e88b04d-44ac-4727-891b-0922a63a25e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "with afw_display.Buffering():\n",
    "    for s in sources:\n",
    "        x = s.getX()\n",
    "        y = s.getY()\n",
    "        if (x > 1500) & (x < 1900) & (y > 1900) & (y < 2300):\n",
    "            afw_display.dot('o', x, y, size=20, ctype='orange')\n",
    "        del x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d232966-e4bd-4046-95ee-e590a2058014",
   "metadata": {},
   "source": [
    "> <img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/dp02_05_fig2.png\" width=300>\n",
    "\n",
    "> Figure 2: A screenshot of what the Firefly interface should look like with the cutout displayed and detected sources overplotted with orange circles.\n",
    "\n",
    "Clean up. Keep `sources` and `result` for the next sections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "574a9b06-cd7e-4ca5-9b34-fe10fca72809",
   "metadata": {},
   "outputs": [],
   "source": [
    "del bbox, cutout, tab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0695c973-e7ce-4c51-9d0e-327a7b0a0847",
   "metadata": {},
   "source": [
    "## 4. Footprints\n",
    "\n",
    "Object footprints are an integral component of the high-level CCD processing tasks (e.g., detection, measurement, and deblending).\n",
    "\n",
    "To quote [Bosch et al. (2017)](https://arxiv.org/pdf/1705.06766.pdf): *\"Footprints record the exact above-threshold detection region on a CCD. These are similar to  SExtractor’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 because they both identify pixels with values above some threshold.\n",
    "\n",
    "The result of the `SourceDetectionTask` stores the footprints associated with detected objects.\n",
    "\n",
    "Grab the above-threshold footprints that were detected, and assign them to the variable `fps`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01a5c82b-9c2b-4f5a-8d0b-0df92182d7a3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fpset = result.positive\n",
    "fps = fpset.getFootprints()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd1400ec-2403-4858-b871-855ebd53723f",
   "metadata": {},
   "source": [
    "Get a rough view of the first source's footprint from its span."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c533489f-d53e-4a72-a4f3-6e0601a20d21",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fps[0].getSpans()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81d27c87-6212-43e3-a965-aabcde5d9921",
   "metadata": {},
   "source": [
    "Above, the footprint is represented by the 1s in the array.\n",
    "\n",
    "Note that the first row of the span array will be the *bottom* row of the image, and so is \"upside down\" compared to the heavy footprint visualization below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ba5221-a8c0-4e67-b803-ee7f09fe154e",
   "metadata": {},
   "source": [
    "### 4.1. Heavy Footprints\n",
    "\n",
    "The footprint span above indicates which pixels are included in the footprint (the 1s).\n",
    "\n",
    "A heavy fooprint includes the pixel value data from the image.\n",
    "\n",
    "Demonstrate that the `fps` footprint is not heavy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5859334b-b7b8-407c-bd77-f26b431aa9cc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fps[0].isHeavy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ba12098-913e-4a07-b77c-79a3794c770b",
   "metadata": {},
   "source": [
    "Make all of the footprints heavy and define `hfps` to hold the heavy footprints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3026cdd5-66e7-4a17-95a1-375f8cb34da7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fpset.makeHeavy(calexp.getMaskedImage())\n",
    "hfps = fpset.getFootprints()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66330675-1a8a-4729-8b14-4b321a7153e1",
   "metadata": {},
   "source": [
    "The `getImageArray` method returns a flattened, 1-dimensional array of pxiels from the footprint.\n",
    "\n",
    "**Option:** Uncomment the line below to print the array and see that it contains pixel values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a69669f9-210c-476e-8ec2-baf9ef3ea3ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# hfps[0].getImageArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5119ea4e-a0aa-470f-9dd7-350eb516d477",
   "metadata": {},
   "source": [
    "Use the span's method to unflatten the image array, and display the heavy footprint ([matplotlib colormap options](https://matplotlib.org/stable/tutorials/colors/colormaps.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6c81694-1701-4f34-b199-dec8808281fc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.imshow(fps[0].getSpans().unflatten(hfps[0].getImageArray()),\n",
    "           cmap='bone', origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62eb82b3-8292-4873-9d56-0673b16644f3",
   "metadata": {},
   "source": [
    "> Figure 3: The heavy footprint of the first source in the `result` array from `SourceDetectionTask`.\n",
    "\n",
    "The heavy footprint also comes with a 1-dimensional mask array.\n",
    "\n",
    "**Option:** Display the 1-d mask array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ad83499-fbb1-446b-96e6-c91543837abe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# hfps[0].getMaskArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ab89f6-8111-410f-a6a3-bf0e4f0a99ed",
   "metadata": {},
   "source": [
    "Review the dictionary for the mask plane, which maps mask value to the reason the pixel is masked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3c49846-d734-445c-895e-4e7700287eaf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "calexp.getMask().getMaskPlaneDict()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9730bce0-bfa6-4b57-b183-b7a21188290c",
   "metadata": {},
   "source": [
    "The values are the exponent of the bitmask.\n",
    "\n",
    "Pixels only marked detected will be 2^5 = 32.\n",
    "\n",
    "Pixels that are both on the edge of the original image and detected will be 2^5 + 2^4 = 48.\n",
    "\n",
    "Visualize mask array.\n",
    "\n",
    "To show the colorbar, create a new axis `cax` on the right side of the image display,\n",
    "with a width of 5% and a padding between display and colorbar of 1%.\n",
    "Add ticks to the colorbar at values of 32 and 48, to correspond to the two\n",
    "mask values of interest (\"detected\" and \"detected + edge\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9b04e62-64da-4f15-afd6-4d162999d02b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "im = plt.imshow(fps[0].getSpans().unflatten(hfps[0].getMaskArray()),\n",
    "                origin='lower')\n",
    "ax = plt.gca()\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"1%\")\n",
    "plt.colorbar(im, cax=cax, ticks=[0, 32, 32+16])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64c4548c-8108-4748-ba28-f7a7b0de2ecb",
   "metadata": {},
   "source": [
    "> Figure 4: Above, the footprint's mask array reveals which pixels are considered (green) but which are also at the edge of the `calexp` (yellow)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81ea7781-6293-4d7a-bd5c-36b412139800",
   "metadata": {},
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "357983f3-3b45-4dad-a5f1-cd19e7c283b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "del fpset, fps, hfps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d884603-978f-4a7f-9988-c9689fa0c7a0",
   "metadata": {},
   "source": [
    "## 5. Forced photometry\n",
    "\n",
    "The LSST Science Pipelines have a task for <a href=\"https://pipelines.lsst.io/modules/lsst.meas.base/tasks_and_algorithms.html#forced-photometry\">forced photometry</a>, called `ForcedMeasurementTask`.\n",
    "\n",
    "This task will use the PSF of the image to measure flux at any user-supplied coordinates.\n",
    "As it is not a *fit*, a flux is returned whether an astronomical source is there or not.\n",
    "Forced photometry can be very useful for including sources below the detection\n",
    "threshold (or non-detections) in analyses.\n",
    "\n",
    "### 5.1. Configure the task\n",
    "\n",
    "To start, create a minimal schema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4bdfdb5-b392-47e1-a5ba-fb9900f0ad94",
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = afwTable.SourceTable.makeMinimalSchema()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d281fae-9f1c-4c30-a798-229dd565af3c",
   "metadata": {},
   "source": [
    "At this time, the `ForcedMeasurementTask` expects a few additional columns related to the centroids\n",
    "and shapes to pre-exist, so add them to the minimal schema.\n",
    "These might change in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8995af12-7ab3-48ea-9cf6-f6702b83e539",
   "metadata": {},
   "outputs": [],
   "source": [
    "alias = schema.getAliasMap() \n",
    "x_key = schema.addField(\"centroid_x\", type=\"D\")\n",
    "y_key = schema.addField(\"centroid_y\", type=\"D\")\n",
    "alias.set(\"slot_Centroid\", \"centroid\")\n",
    "\n",
    "xx_key = schema.addField(\"shape_xx\", type=\"D\")\n",
    "yy_key = schema.addField(\"shape_yy\", type=\"D\")\n",
    "xy_key = schema.addField(\"shape_xy\", type=\"D\")\n",
    "alias.set(\"slot_Shape\", \"shape\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfb37786-17d5-41c1-a2a8-ada0a065ccfe",
   "metadata": {},
   "source": [
    "Add a type flag, which will be 0 for sub-threshold (faint) sources and 1 for random \"off-source\" coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa0c1dc-e259-45a3-b497-05b55ca06663",
   "metadata": {},
   "outputs": [],
   "source": [
    "type_key = schema.addField(\"type_flag\", type=\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ad48ac4-8167-4da2-86f7-9959e8c9a231",
   "metadata": {},
   "source": [
    "Create an empty `forcedSource` catalog using the schema. It will have 0 rows, to start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f7ef27a-c22e-4f32-a022-b710e3d8c5c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "forcedSource = afwTable.SourceCatalog(schema)\n",
    "print(len(forcedSource))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1911ae25-1e09-4f05-a642-859015965ed1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-17T19:08:40.605507Z",
     "iopub.status.busy": "2024-05-17T19:08:40.604853Z",
     "iopub.status.idle": "2024-05-17T19:08:40.817706Z",
     "shell.execute_reply": "2024-05-17T19:08:40.816743Z",
     "shell.execute_reply.started": "2024-05-17T19:08:40.605484Z"
    }
   },
   "source": [
    "Configure the `ForcedMeasurementTask`. At this time, few plugins beyond the default are required, as is the `doReplaceWithNoise` configuration. These might change in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62dceb8b-6be0-4e45-bd36-7b98324b028a",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = ForcedMeasurementTask.ConfigClass()\n",
    "config.copyColumns = {}\n",
    "config.plugins.names = [\n",
    "    \"base_TransformedCentroid\",\n",
    "    \"base_PsfFlux\",\n",
    "    \"base_TransformedShape\"\n",
    "]\n",
    "config.doReplaceWithNoise = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9571486-7023-47f0-a937-59c357580608",
   "metadata": {},
   "source": [
    "Define the task using the schema and configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e31f102-87dd-434f-bdd7-089d338b4280",
   "metadata": {},
   "outputs": [],
   "source": [
    "forcedMeasurementTask = ForcedMeasurementTask(schema, config=config)\n",
    "del config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4901d6a-68f7-4edc-9bb5-c21005561819",
   "metadata": {},
   "source": [
    "### 5.2. Define the coordinates\n",
    "\n",
    "Scientific motivations for forced photometry include measurement of sub-threshold (faint) sources or a \n",
    "need for an upper-limit on the flux of a time-varying source.\n",
    "\n",
    "The potential drawback of forced photometry is that all the flux within the area of the PSF is included,\n",
    "which means nearby sources can contaminate (add excess flux to) the forced photometry.\n",
    "This will be demonstrated in Section 5.4.\n",
    "\n",
    "For this example, create a list of sub-threshold (faint) coordinates and random \"off-source\" coordinates which are not near any detected sources.\n",
    "Use a region within 2 arcmin of the center of the field.\n",
    "\n",
    "#### 5.2.1. Sub-threshold (faint) sources\n",
    "\n",
    "Create `numpy` arras of the right ascension and declination, x- and y-pixel coordinates, and instrumental flux and its error\n",
    "from the list of detected sources. Calculate the signal-to-noise ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07645563-4f0a-4950-9ff1-c3cbaf115aac",
   "metadata": {},
   "outputs": [],
   "source": [
    "srcs_ra = np.asarray(np.rad2deg(sources['coord_ra']), dtype='float')\n",
    "srcs_dec = np.asarray(np.rad2deg(sources['coord_dec']), dtype='float')\n",
    "srcs_xpix = np.asarray(sources['base_SdssCentroid_x'], dtype='float')\n",
    "srcs_ypix = np.asarray(sources['base_SdssCentroid_y'], dtype='float')\n",
    "srcs_flux = np.asarray(sources['base_PsfFlux_instFlux'], dtype='float')\n",
    "srcs_fluxe = np.asarray(sources['base_PsfFlux_instFluxErr'], dtype='float')\n",
    "srcs_snr = srcs_flux/srcs_fluxe"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d308f9a-3723-42fa-b8b7-2e410208f8fd",
   "metadata": {},
   "source": [
    "Use the `bbox` and `wcs` to get the central coordinates of the `calexp`, convert from pixels to sky coordinates,\n",
    "and calculate offsets (`s_dra` and `s_ddec`) between the sources and the image center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4992ebe-3db0-4247-b9ba-3211d68acfd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "bbox = calexp.getBBox()\n",
    "wcs = calexp.getWcs()\n",
    "center = wcs.pixelToSky(bbox.centerX, bbox.centerX)\n",
    "cra = center.getRa().asDegrees()\n",
    "cdec = center.getDec().asDegrees()\n",
    "srcs_dra = np.abs(srcs_ra - cra)\n",
    "srcs_ddec = np.abs(srcs_dec - cdec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "414b622b-5881-49ae-83be-2a63b7b616a6",
   "metadata": {},
   "source": [
    "Define the maximum offset from the image center to be 2 arcminutes (in units of degrees).\n",
    "\n",
    "Find all sources with 4 $<$ SNR $<$ 5 (i.e., sub-threshold sources) within the central region.\n",
    "\n",
    "Define new arrays to hold only the information of the sources to be used for forced photometry, and clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31cf3c7b-96e5-4313-a375-b703f4de992d",
   "metadata": {},
   "outputs": [],
   "source": [
    "dmax = 2.0 / 60.0\n",
    "tx = np.where((srcs_snr > 4) & (srcs_snr < 5)\n",
    "              & (srcs_dra < dmax) & (srcs_ddec < dmax))[0]\n",
    "print(len(tx))\n",
    "\n",
    "s_ra = srcs_ra[tx]\n",
    "s_dec = srcs_dec[tx]\n",
    "s_xpix = srcs_xpix[tx]\n",
    "s_ypix = srcs_ypix[tx]\n",
    "s_flux = srcs_flux[tx]\n",
    "\n",
    "del tx, dmax\n",
    "del srcs_ra, srcs_dec, srcs_xpix, srcs_ypix, srcs_flux, srcs_fluxe, srcs_snr\n",
    "del srcs_dra, srcs_ddec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2821f9e9-f511-4136-aca8-9078c94b2eb4",
   "metadata": {},
   "source": [
    "Add the coordinates of these 32 sub-threshold sources to the `forcedSource` catalog.\n",
    "\n",
    "Sky coordinates must be in radians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da08e895-5dc4-4704-81e0-f34183fbe329",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(s_ra)):\n",
    "    sourceRec = forcedSource.addNew()\n",
    "    coord = geom.SpherePoint(np.radians(s_ra[i]) * geom.radians,\n",
    "                             np.radians(s_dec[i]) * geom.radians)\n",
    "    sourceRec.setCoord(coord)\n",
    "    sourceRec[x_key] = s_xpix[i]\n",
    "    sourceRec[y_key] = s_ypix[i]\n",
    "    sourceRec[type_key] = 0\n",
    "    del sourceRec, coord\n",
    "\n",
    "print(len(forcedSource))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7f21859-39f5-4298-bfdb-c17cd59b58e1",
   "metadata": {},
   "source": [
    "#### 5.2.2. Random \"off-source\" coordinates\n",
    "\n",
    "Create 32 random coordinates in the central region of the image. Ensure they are offset by at least 40 pixels from the nearest detected source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f09aa4-6b59-42ef-958c-36c2ff64d5cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "pixelScale = calexp.getWcs().getPixelScale().asArcseconds()\n",
    "radius_pix = int(2.0 * 60.0 / pixelScale)\n",
    "center_x = bbox.centerX\n",
    "center_y = bbox.centerY\n",
    "\n",
    "tempx = []\n",
    "tempy = []\n",
    "while len(tempx) < 32:\n",
    "    r = np.random.choice(radius_pix)\n",
    "    theta = np.random.rand() * 2 * np.pi\n",
    "    x = center_x + r * np.cos(theta)\n",
    "    y = center_y + r * np.sin(theta)\n",
    "    off = np.min(np.sqrt((x - s_xpix)**2 + (y - s_ypix)**2))\n",
    "    if off > 40:\n",
    "        tempx.append(x)\n",
    "        tempy.append(y)\n",
    "    del r, theta, x, y, off\n",
    "\n",
    "r_xpix = np.asarray(tempx)\n",
    "r_ypix = np.asarray(tempy)\n",
    "del pixelScale, radius_pix, center_x, center_y, tempx, tempy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19897e58-7209-469b-899a-d03f93b598f0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-02T18:17:35.326935Z",
     "iopub.status.busy": "2024-06-02T18:17:35.326375Z",
     "iopub.status.idle": "2024-06-02T18:17:35.329765Z",
     "shell.execute_reply": "2024-06-02T18:17:35.329411Z",
     "shell.execute_reply.started": "2024-06-02T18:17:35.326917Z"
    }
   },
   "source": [
    "Use the image's WCS to convert the random pixels to sky coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "994db175-7612-4ed2-b2f5-33c69a6db186",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_ra = np.zeros(len(r_xpix), dtype='float')\n",
    "r_dec = np.zeros(len(r_ypix), dtype='float')\n",
    "for i in range(len(r_xpix)):\n",
    "    r_skycoord = wcs.pixelToSky(r_xpix[i], r_ypix[i])\n",
    "    r_ra[i] = r_skycoord.getRa().asDegrees()\n",
    "    r_dec[i] = r_skycoord.getDec().asDegrees()\n",
    "    del r_skycoord"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "861596e1-05bf-4712-9166-e46caa4f021b",
   "metadata": {},
   "source": [
    "Add the random \"off-source\" coordinates to the `forcedSource` catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd146247-392c-49a1-a6e3-7d95318871d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(r_xpix)):\n",
    "    sourceRec = forcedSource.addNew()\n",
    "    coord = geom.SpherePoint(np.radians(r_ra[i]) * geom.radians,\n",
    "                             np.radians(r_dec[i]) * geom.radians)\n",
    "    sourceRec.setCoord(coord)\n",
    "    sourceRec[x_key] = r_xpix[i]\n",
    "    sourceRec[y_key] = r_ypix[i]\n",
    "    sourceRec[type_key] = 1\n",
    "    del sourceRec, coord\n",
    "\n",
    "print(len(forcedSource))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64ea2ae8-2f22-4d41-8373-7d9cc927b92d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-17T19:08:17.631865Z",
     "iopub.status.busy": "2024-05-17T19:08:17.630985Z",
     "iopub.status.idle": "2024-05-17T19:08:18.126891Z",
     "shell.execute_reply": "2024-05-17T19:08:18.125931Z",
     "shell.execute_reply.started": "2024-05-17T19:08:17.631830Z"
    }
   },
   "source": [
    "**Option:** View the `forcedSource` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2a62d43-6237-40ed-8744-d227ab59e74e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# forcedSource"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e56a38b-012c-4843-be40-68b6f939d0cb",
   "metadata": {},
   "source": [
    "### 5.3. Run the task\n",
    "\n",
    "Create an empty catalog to store forced measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "200df4af-dbb1-4863-8f5f-370fe7446f6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "forcedMeasCat = forcedMeasurementTask.generateMeasCat(calexp, forcedSource,\n",
    "                                                      calexp.getWcs())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8948d9d0-b9dc-49f6-bc2c-eda7e03b62ac",
   "metadata": {},
   "source": [
    "Option to review the `forcedMeasCat`.\n",
    "Before the task is run, all the columns are empty (or default values)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fe1e6fd-2de8-42d9-88dc-17aa101ab326",
   "metadata": {},
   "outputs": [],
   "source": [
    "# forcedMeasCat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c769ba8a-0a4e-45af-af0b-2bb59db65a8e",
   "metadata": {},
   "source": [
    "Run the forced measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09495337-4887-49a3-ac06-cbcf4a6bf46b",
   "metadata": {},
   "outputs": [],
   "source": [
    "forcedMeasurementTask.run(forcedMeasCat, calexp, forcedSource, calexp.getWcs())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "031fea29-2c33-4c9a-8ed2-565c0c33b558",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-17T19:10:41.842379Z",
     "iopub.status.busy": "2024-05-17T19:10:41.841667Z",
     "iopub.status.idle": "2024-05-17T19:10:42.055329Z",
     "shell.execute_reply": "2024-05-17T19:10:42.054237Z",
     "shell.execute_reply.started": "2024-05-17T19:10:41.842355Z"
    }
   },
   "source": [
    "### 5.4. Explore the results\n",
    "\n",
    "Get the results as an `astropy` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8af77f19-8d84-4d16-81ea-6cfdccf5bce3",
   "metadata": {},
   "outputs": [],
   "source": [
    "table = forcedMeasCat.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77a511f2-e3a5-497a-a38f-0708d161c948",
   "metadata": {},
   "source": [
    "Add the `type_flag` column from `forcedSource` to `table`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65bf0bd5-63d8-408c-bc93-45db98abf9d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = forcedSource.asAstropy()\n",
    "table['type_flag'] = temp['type_flag']\n",
    "del temp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfd34e8d-28f1-4a02-a63d-dd6bf12e0adf",
   "metadata": {},
   "source": [
    "**Option:** View the filled `forcedMeasCat` results as the `astropy` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0db517b7-b4e3-4206-b63b-ec0622c8dae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fee6cfc5-ed4a-4a8e-ba1b-566a39c451c7",
   "metadata": {},
   "source": [
    "Display the full `calexp` in frame 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eb05112-4222-405b-9746-5f2485678b2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display = afwDisplay.Display(frame=1)\n",
    "afw_display.mtv(calexp)\n",
    "afw_display.setMaskTransparency(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0427a846-b1d8-4c5c-8a6c-d5269949160f",
   "metadata": {},
   "source": [
    "Overplot the coordinates used for forced photometry on the image.\n",
    "\n",
    "Use circles to mark sub-threshold (faint) source coodinates (orange) and random \"off-source\" coordinates (cyan), showing the locations where forced photometry was done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4b7ae7b-040f-435d-8f96-f6ee484d05a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "with afw_display.Buffering():\n",
    "    for i in range(len(table)):\n",
    "        clr = 'orange'\n",
    "        if table['type_flag'][i] == 1:\n",
    "            clr = 'cyan'\n",
    "        afw_display.dot('o', table['slot_Centroid_x'][i],\n",
    "                        table['slot_Centroid_y'][i],\n",
    "                        size=20, ctype=clr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcd84e36-d409-48aa-b20a-a48d6a9aca75",
   "metadata": {},
   "source": [
    "Zoom-in on the displayed image in Firefly to further explore the\n",
    "locations where forced photometry was performed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "083760c2-3e95-4a09-a3db-1d7839d78c53",
   "metadata": {},
   "source": [
    "Compare the forced photometry fluxes with the measured fluxes for sub-threshold (faint) sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d0378dc-fd9c-4540-a8f7-bd832606b096",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "plt.plot([2500, 3000], [2500, 3000], lw=0.5, color='black')\n",
    "for i in range(32):\n",
    "    plt.plot(s_flux[i], table['base_PsfFlux_instFlux'][i],\n",
    "             'o', color='darkorange', mew=0, alpha=0.5)\n",
    "plt.xlabel('detection flux')\n",
    "plt.ylabel('forced PSF flux')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83feb3cc-85f0-44f4-803b-93efb1b4ef17",
   "metadata": {},
   "source": [
    "> Figure 5: Above, the original detection instrumental flux (in counts) versus the forced PSF flux measurements for the 32 sources which were originally detected with 4 $<$ SNR $<$ 5. For most sources the fluxes agree, but in some cases the forced PSF flux is greater.\n",
    "\n",
    "Cases in which the forced PSF flux is greater is likely due to contamination from a nearby source. This contamination is mitigated when source deblending and measurement are done together, but forced photometry includes all the flux in the area of the PSF.\n",
    "\n",
    "Again overplot sources on the image, but color as yellow the ones with a flux excess greater than 20 counts. Notice that for at least five, the contaminating source is visible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a036eff-f481-4383-b7b3-734e8b9f87f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "with afw_display.Buffering():\n",
    "    for i in range(len(table)):\n",
    "        if table['type_flag'][i] == 0:\n",
    "            if table['base_PsfFlux_instFlux'][i] > s_flux[i] + 20:\n",
    "                clr = 'yellow'\n",
    "                afw_display.dot('o', table['slot_Centroid_x'][i],\n",
    "                                table['slot_Centroid_y'][i],\n",
    "                                size=20, ctype=clr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "080737b8-a576-4109-8816-fa4f79ac754c",
   "metadata": {},
   "source": [
    "> <img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/dp02_05_fig6.png\" width=300>\n",
    "\n",
    "> Figure 6: A screenshot of the calexp displayed with the locations of forced photometry marked by circles. Coordinates of sub-threshold detections are marked with orange, or yellow if their forced PSF flux was significantly in excess of the deblended flux. Random \"off-source\" coordinates are marked in cyan.\n",
    "\n",
    "As a final exploration, plot the histograms of the forced PSF fluxes for the sub-threshold (faint) sources and\n",
    "the random \"off-source\" coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7956122-9d96-4421-9e0f-9bab0aedf369",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "tx1 = np.where(table['type_flag'] == 0)[0]\n",
    "tx2 = np.where(table['type_flag'] == 1)[0]\n",
    "plt.hist(table['base_PsfFlux_instFlux'][tx1], bins=10, histtype='step',\n",
    "         color='darkorange', label='sub-threshold source')\n",
    "plt.hist(table['base_PsfFlux_instFlux'][tx2], bins=40, histtype='step',\n",
    "         color='dodgerblue', label='random coordinate')\n",
    "plt.xlabel('forced PSF flux')\n",
    "plt.ylabel('number of forced sources')\n",
    "plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e5d0e06-4abf-4c52-9226-52c7a717a28e",
   "metadata": {},
   "source": [
    "> Figure 7: Above, histograms of the forced PSF flux for random \"off-source\" coordinates (blue) and for the coordinates of detected sub-threshold sources (orange). As expected, the \"off-source\" fluxes are centered on zero flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21f092f2-8804-4c0e-8942-da7b2c53b0f3",
   "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
}