2017_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the LSST DM Stack in Python\n",
    "\n",
    "This tutorial focuses on using the DM stack in Python.\n",
    "The tasks we will go through in this notebook are more commonly executed as part of larger processing pipelines.\n",
    "We execute them in a notebook here to give a peek into the pipelines and hopefully demonstrate how these pieces of the pipelines can be configured.\n",
    "\n",
    "More notebook examples can be found here: https://github.com/RobertLuptonTheGood/notebooks/tree/master/Demos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Data Repository\n",
    "\n",
    "Instead of operating directly on files and directories, we interact with on-disk data products via an abstraction layer called the *data butler*.\n",
    "The butler operates on *data repositories*.\n",
    "There is a small data repository generated from HyperSuprimeCam (HSC) public data located at `/project/shared/data/RSP_CHECK_REPO/` on all deployments of the Rubin Science Platform.\n",
    "In this tutorial we'll focus on processing an individual image, and that's all this particular subset will support."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Access with Butler\n",
    "\n",
    "Data repositories are managed by an object called `Butler`.\n",
    "To retrieve a dataset from the `Butler`, we start by constructing one pointing to the repository containing the data used in this tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.daf.butler import Butler\n",
    "butler = Butler('/project/shared/data/RSP_CHECK_REPO/butler.yaml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then call `get` with the name and data ID of the dataset.\n",
    "The name of the image that's saved directly after instrument signature removal (ISR) is `postISRCCD`.\n",
    "We will put this and the code necessary to extract a sub region in a function that we can call later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regen_exposure():\n",
    "    exposure = butler.get(\"postISRCCD\", exposure=1228, detector=87, instrument='HSC',\n",
    "                          collections='HSC/defaults')\n",
    "\n",
    "    bbox = exposure.getBBox()\n",
    "    bbox.grow(-bbox.getDimensions()//3)  # box containing the central third (in each dimension)\n",
    "    bbox.grow(-Extent2I(0, 400))  # make it a bit smaller in x\n",
    "    # exposure[bbox] would also work here because exposure.getXY0() == (0, 0),\n",
    "    # but it's dangerous in general because it ignores that origin.\n",
    "    sub = Exposure(exposure, bbox=bbox, dtype=exposure.dtype, deep=False)\n",
    "    return exposure, sub, bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image, Boxes, and (Crude) Image Display\n",
    "\n",
    "A full 2k x 4k HSC CCD is a pretty big image to display when you don't have specialized display code.  The DM stack does have specialized display code, but it either requires DS9 (which requires some ssh tunnels to use with data living on a server) or a Firefly server installation. For this tutorial, we'll just throw together a naive matplotlib display function, and create a view to a subimage that we'll display instead of the full image.\n",
    "\n",
    "This section features a few of our most important class objects:\n",
    "\n",
    "- `lsst.afw.image.Exposure` is an image object that actually holds three image planes: the science image (`Exposure.image`), an image of variance in every pixel (`Exposure.variance`), an integer bit mask (`Exposure.mask`).  It also holds a lot of more complex objects that characterize the image, such as a point-spread function (`lsst.afw.detection.Psf`) and world-coordinate system (`lsst.afw.image.Wcs`).  Most of these objects aren't filled in yet, because all we've run so far is ISR.  It doesn't generally make sense to perform mathematical operations (i.e. addition) on `Exposure`s, because those operations aren't always well-defined on the more complex objects.  You can get a `MaskedImage` object with the same image, mask, and variance planes that does support mathematical operations but doesn't contain `Psf`s and `Wcs`s (etc) with `Exposure.maskedImage`.\n",
    "\n",
    "- The `Exposure.image` and `Exposure.variance` properties return `lsst.afw.image.Image` objects.  These have a `.array` property that returns a `numpy.ndarray` view to the `Image`'s pixels.  Conceptually, you should think of an `Image` as just a `numpy.ndarray` with a possibly nonzero origin.\n",
    "\n",
    "- The `Exposure.mask` property returns a `lsst.afw.image.Mask` object, which behaves like an `Image` with a dictionary-like object that relates string labels to bit numbers.\n",
    "\n",
    "- All of these image-like objects have a `getBBox()` method, which returns a `lsst.afw.geom.Box2I`.  The minimum and maximum points of a `Box2I` are specified in integers that correspond to the *centers* of the lower-left and upper-right pixels in the box, but the box conceptually contains the entirety of those pixels.  To get a box with a floating-point representation of the same boundary for the `extent` argument to `imshow` below, we construct a `Box2D` from the `Box2I`.\n",
    "\n",
    "- `Point2I` and `Extent2I` are used to represent absolute positions and offsets between positions as integers (respectively).  These have floating-point counterparts `Point2D` and `Extent2D`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.geom import Box2D, Box2I, Point2I, Extent2I\n",
    "from lsst.afw.image import Exposure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams[\"figure.figsize\"] = (8, 6)\n",
    "matplotlib.rcParams[\"font.size\"] = 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def display(image, mask=None, colors=None, alpha=0.40, **kwds):\n",
    "    box = Box2D(image.getBBox())\n",
    "    extent = (box.getMinX(), box.getMaxX(), box.getMinY(), box.getMaxY())\n",
    "    kwds.setdefault(\"extent\", extent)\n",
    "    kwds.setdefault(\"origin\", \"lower\")\n",
    "    kwds.setdefault(\"interpolation\", \"nearest\")\n",
    "    matplotlib.pyplot.imshow(image.array, **kwds)\n",
    "    kwds.pop(\"vmin\", None)\n",
    "    kwds.pop(\"vmax\", None)\n",
    "    kwds.pop(\"norm\", None)\n",
    "    kwds.pop(\"cmap\", None)\n",
    "    if mask is not None:\n",
    "        for plane, color in colors.items():\n",
    "            array = np.zeros(mask.array.shape + (4,), dtype=float)\n",
    "            rgba = np.array(matplotlib.colors.hex2color(matplotlib.colors.cnames[color]) + (alpha, ),\n",
    "                            dtype=float)\n",
    "            np.multiply.outer((mask.array & mask.getPlaneBitMask(plane)).astype(bool), rgba, out=array)\n",
    "            matplotlib.pyplot.imshow(array, **kwds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now here's the (cutout) of the detrended image.  I've cheated in setting the scale by looking at the background level in advance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(sub.image, vmin=3300, vmax=3500, cmap=matplotlib.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background Subtraction and Task Configuration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step we usually take is to estimate and subtract the background, using `lsst.meas.algorithms.SubtractBackgroundTask`.  This is a regular `Task`, not a `CmdLineTask`, and hence we'll just pass it our `Exposure` object (it operates in-place) instead of a `Butler`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SubtractBackgroundTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgConfig = SubtractBackgroundTask.ConfigClass()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to get fun & terrible results!\n",
    "bkgConfig.useApprox = False\n",
    "bkgConfig.binSize = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pattern for configuration here is the same as it was for `SubaruIsrTask`, but here we're setting values directly instead of loading a configuration file from the `obs_subaru` camera-specialization package.  The `config` object here is an instance of a class that inherits from `lsst.pex.config.Config` that contains a set of `lsst.pex.config.Field` objects that define the options that can be modified.  Each `Field` behaves more or less like a Python `property`, and you can get information on all of the fields in a config object by either using `help`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(bkgConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SubtractBackgroundTask.ConfigClass.algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgTask = SubtractBackgroundTask(config=bkgConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgResult = bkgTask.run(exposure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(sub.image, vmin=-0.5, vmax=100, cmap=matplotlib.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've run through all of these steps after executing the cell that warns about terrible results, you should notice that the galaxy in the upper right has been oversubtracted.\n",
    "\n",
    "**EXERCISE**: Before continuing on, re-load the exposure from disk, reset the configuration and `Task` instances, and re-run without executing the cell that applies bad values to the config, all by just re-executing the right cells above.  You should end up an image in which the upper-right galaxy looks essentially the same as it does in the image before we subtracted the background."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ANSWER:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could construct a new instance of `SubtractBackgroundTask.ConfigClass` and pass it unmodified to the constructor (and that's what you'd get by strictly copying existing cells), but it also works to just construct the task with no arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgTask = SubtractBackgroundTask()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgResult = bkgTask.run(exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installing an Initial-Guess PSF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most later processing steps require a PSF model, which is represented by a `Psf` object that's attached to the `Exposure`.  For now, we'll just make a Gaussian PSF with some guess at the seeing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SingleGaussianPsf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FWHM_TO_SIGMA = 1.0/(2*np.sqrt(2*np.log(2)))\n",
    "PIXEL_SCALE = 0.168  # arcsec/pixel\n",
    "SEEING = 0.7         # FWHM in arcsec\n",
    "sigma = FWHM_TO_SIGMA*SEEING/PIXEL_SCALE\n",
    "width = int(sigma*3)*2 + 1\n",
    "psf = SingleGaussianPsf(width, width, sigma=sigma)\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `Psf` object can basically just do one thing: it can return an image of itself at a point.  `SingleGaussianPsf` represents a constant PSF, so it always returns the same image, regardless of the point you give it.\n",
    "\n",
    "But there are two ways to evaluate a `Psf` at a point.  If you want an image centered on the middle pixel, and that middle pixel to be the origin - what you'd usually want if you're going to convolve the PSF with another model - use `computeKernelImage(point)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.geom import Point2D\n",
    "display(psf.computeKernelImage(Point2D(60.5, 7.2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to compare the PSF to a star at the exact same position, use `computeImage(point)`.  That will shift the image returned by `computeKernelImage(point)` to the right sub-pixel offset, and update the origin of the image to take care of the rest, so you end up with a postage stamp in the same coordinate system as the original image where the star is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(psf.computeImage(Point2D(60.5, 7.2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Removing Cosmic Rays\n",
    "\n",
    "Cosmic rays are detected and interpolated by `RepairTask`, which also sets mask planes to indicate where the cosmic rays were (\"CR\") and which pixels were interpolated (\"INTERP\"; this may happen due to saturation or bad pixels as well).  Because we're just using the default configuration, we can skip creating a config object and just construct the `Task` with no arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.pipe.tasks.repair import RepairTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repairTask = RepairTask()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repairTask.run(exposure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(sub.image, mask=sub.mask, colors={\"CR\": \"red\"},\n",
    "        vmin=-0.5, vmax=100, alpha=0.8, cmap=matplotlib.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detecting Sources\n",
    "\n",
    "Unlike the other `Task`s we've dealt with so far, `SourceDetectionTask` creates a `SourceCatalog` in addition to updating the image (all it does to the image is add a \"DETECTED\" mask plane).  All `Task`s that work with catalogs need to be initialized with a `lsst.afw.table.Schema` object, to which the `Task` will add the fields necessary to store its outputs.  A `SourceCatalog`'s `Schema` cannot be modified after the `SourceCatalog` has been constructed, which means it's necessary to construct all `Schema`-using `Task`s before actually running any of them.\n",
    "\n",
    "Each record in the catalog returned by `SourceDetectionTask` has a `Footprint` object attached to it.  A `Footprint` represents the approximate region covered by a source in a run-length encoding data structure.  It also contains a list of peaks found within that region.  The \"DETECTED\" mask plane is set to exactly the pixels covered by any `Footprint` in the returned catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SourceDetectionTask\n",
    "from lsst.afw.table import SourceTable, SourceCatalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = SourceTable.makeMinimalSchema()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detectTask = SourceDetectionTask(schema=schema)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A SourceTable is really just a factory object for records; don't confuse it with SourceCatalog, which is\n",
    "# usually what you want.  But a SourceTable *is* what SourceDetectionTask wants here.\n",
    "table = SourceTable.make(schema)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detectResult = detectTask.run(table, exposure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(sub.image, mask=sub.mask, colors={\"DETECTED\": \"blue\"}, vmin=-0.5, vmax=100, cmap=matplotlib.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deblending\n",
    "\n",
    "Deblending attempts to separate detections with multiple peaks into separate objects.  We keep all of the original sources in the `SourceCatalog` (called `parent`s) when we deblend, but for each `parent` source that contains more than one peak, we create a new record (called a `child`) for each of those peaks.  The `Footprint`s attached to the `child` objects are instances of a subclass called `HeavyFootprint`, which include new deblended pixel values as well as the region description.  These can be used by calling `insert` to replace an `Image`'s pixels with the `HeavyFootprint`'s pixels.\n",
    "\n",
    "**EXERCISE**: This section will not run if the cells are executed naively in order.  At some point you'll have to go re-execute one or more cells in the previous section to get the right behavior.  Which one(s)?  Why?  Copy those cells here (in the right places) when you figure it out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ANSWER:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.meas.deblender import SourceDeblendTask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It isn't necessary to re-load the `Exposure` from disk (and re-do background subtraction) to get this exercise to run, but not doing so will result in re-subtracting the residual background repeatedly every time `SourceDetectionTask` is run, which means the results will change slightly each time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It isn't necessary to construct all of these tasks before running any of them (just the ones that take and modify a schema), but it's easier to get things right if you follow that principle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = SourceTable.makeMinimalSchema()\n",
    "bkgTask = SubtractBackgroundTask()\n",
    "repairTask = RepairTask()\n",
    "detectTask = SourceDetectionTask(schema=schema)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deblendTask = SourceDeblendTask(schema=schema)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run the tasks from the previous sections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgResult = bkgTask.run(exposure)\n",
    "repairTask.run(exposure)\n",
    "table = SourceTable.make(schema)\n",
    "detectResult = detectTask.run(table, exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, finally, we can run the task introduced in this section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = detectResult.sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deblendTask.run(exposure, catalog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To inspect some deblender outputs, we'll start by finding some parent objects that were deblended into multiple children, by looking at the `deblend_nChild` field (which was added to the `Schema` when we constructed the `SourceDeblendTask`, and populated when we called `run`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find some blended sources inside the subimage:\n",
    "blendParents = []\n",
    "for record in catalog:\n",
    "    if record.get(\"deblend_nChild\") > 0 and bbox.contains(record.getFootprint().getBBox()):\n",
    "        blendParents.append(record)\n",
    "# Sort by peak brightness so we can look at something with decent S/N\n",
    "blendParents.sort(key=lambda r: -r.getFootprint().getPeaks()[0].getPeakValue())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.image import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image of the parent object is just the original image, but we'll cut out just the region inside its `Footprint`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blendParentImage = Image(exposure.image, bbox=blendParents[0].getFootprint().getBBox(),\n",
    "                         deep=True, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll insert the deblended child pixels into blank images of the same size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blendChildImages = []\n",
    "for blendChild in catalog.getChildren(blendParents[0].getId()):\n",
    "    image = Image(blendParentImage.getBBox(), dtype=np.float32)\n",
    "    blendChild.getFootprint().insert(image)\n",
    "    blendChildImages.append(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nSubPlots = len(blendChildImages) + 1\n",
    "nCols = 3\n",
    "nRows = nSubPlots//nCols + 1\n",
    "matplotlib.pyplot.subplot(nRows, nCols, 1)\n",
    "display(blendParentImage, vmin=-0.5, vmax=100, cmap=matplotlib.cm.gray)\n",
    "for n, image in enumerate(blendChildImages):\n",
    "    matplotlib.pyplot.subplot(nRows, nCols, n + 2)\n",
    "    display(image, vmin=-0.5, vmax=100, cmap=matplotlib.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurement\n",
    "\n",
    "`SingleFrameMeasurementTask` is typically responsible for adding most fields to a `SourceCatalog`.  It runs a series of plugins that make different measurements (you can configure them with the `.plugins` dictionary-like field on its config object, and control which are run with `.names`).  If the deblender has been run first, it will measure child objects using their deblended pixels.\n",
    "\n",
    "**EXERCISE**: Like the Deblending section, you'll have to re-execute some previous cells somewhere in this section to get the right behavior.  Copy those cells into the right places here once you've gotten it working."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ANSWER:** As in the Deblending section exercise, we reload the exposure, create all of the tasks, and then run all of the tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.meas.base import SingleFrameMeasurementTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = SourceTable.makeMinimalSchema()\n",
    "bkgTask = SubtractBackgroundTask()\n",
    "repairTask = RepairTask()\n",
    "detectTask = SourceDetectionTask(schema=schema)\n",
    "deblendTask = SourceDeblendTask(schema=schema)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measureConfig = SingleFrameMeasurementTask.ConfigClass()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What measurements are configured to run\n",
    "print(measureConfig.plugins.names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import an extension module that adds a new measurement\n",
    "import lsst.meas.extensions.photometryKron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What measurements *could* be configured to run\n",
    "print(list(measureConfig.plugins.keys()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure the new measurement to run\n",
    "measureConfig.plugins.names.add(\"ext_photometryKron_KronFlux\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measureTask = SingleFrameMeasurementTask(schema=schema, config=measureConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgResult = bkgTask.run(exposure)\n",
    "repairTask.run(exposure)\n",
    "table = SourceTable.make(schema)\n",
    "detectResult = detectTask.run(table, exposure)\n",
    "catalog = detectResult.sources\n",
    "deblendTask.run(exposure, catalog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measureTask.run(catalog, exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll show some of the results of measurement by overlaying the measured ellipses on the image.\n",
    "\n",
    "The shapes and centroids we use here (by calling `record.getX()`, `record.getY()`, `record.getShape()`) are aliases  (called \"slots\") to fields with longer names that are our recommended measurements for these quantities.  You can see the set of aliases by printing the schema (see next section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.geom.ellipses import Axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(sub.image, mask=sub.mask, colors={\"DETECTED\": \"blue\"}, vmin=-0.5, vmax=100, cmap=matplotlib.cm.gray)\n",
    "\n",
    "for record in catalog:\n",
    "    if record.get(\"deblend_nChild\") != 0:\n",
    "        continue\n",
    "    axes = Axes(record.getShape())   # convert to A, B, THETA parameterization\n",
    "    axes.scale(2.0)  # matplotlib uses diameters, not radii\n",
    "    # Don't generate ellipses for records with nans\n",
    "    nans = np.isnan([record.getX(), record.getY(), axes.getA(), axes.getB(), axes.getTheta()])\n",
    "    if sum(nans) > 0:\n",
    "        continue\n",
    "    patch = matplotlib.patches.Ellipse(\n",
    "               (record.getX(), record.getY()),\n",
    "               axes.getA(), axes.getB(),\n",
    "               angle=axes.getTheta() * 180.0 / np.pi,\n",
    "               fill=False, edgecolor=\"green\"\n",
    "            )\n",
    "    matplotlib.pyplot.gca().add_patch(patch)\n",
    "matplotlib.pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working With Catalogs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the schema:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(catalog.getSchema())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get arrays of columns (requires the catalog to be continguous in memory, which we can guarantee with a deep copy):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = catalog.copy(deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psfFlux = catalog[\"base_PsfFlux_instFlux\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that boolean values are stored in `Flag` columns, which are packed into bits.  Unlike other column types, when you get an array of a `Flag` column, you get a copy, not a view."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `Key` objects instead of strings to do fast repeated access to fields when iterating over records:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = catalog.getSchema().find(\"deblend_nChild\").key\n",
    "deblended = [record for record in catalog if record.get(key) == 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also get `dict` version of a subset of a `Schema`, a `Catalog`, or a `Record` by calling either `extract` methods with a glob:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog[0].extract(\"base_PsfFlux_*\")  # or regex='...'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For `Record`s, the dict values are just the values of the fields, and for `Catalogs`, they're `numpy.ndarray` columns.  For `Schema`s they're `SchemaItem`s, which behave liked a named tuple containing a `Key` and a `Field`, which contains more descriptive information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get an Astropy view of the catalog (from which you can make a Pandas view):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = catalog.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can find some reference documentation for the catalog library [here](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/afw_table.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "1) Make a scatter plot of Kron Flux vs. Psf Flux."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's not a pretty plot, with so few points and no calibration, but it's pretty simple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.pyplot.scatter(catalog[\"slot_PsfFlux_instFlux\"], catalog[\"ext_photometryKron_KronFlux_instFlux\"], alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Write a single function that performs all of the above steps on a post-ISR `Exposure` object, modifying the `Exposure` in-place and returning a new `SourceCatalog` with a complete set of measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def processExposure1(exposure):\n",
    "    schema = SourceTable.makeMinimalSchema()\n",
    "    bkgTask = SubtractBackgroundTask()\n",
    "    repairTask = RepairTask()\n",
    "    detectTask = SourceDetectionTask(schema=schema)\n",
    "    deblendTask = SourceDeblendTask(schema=schema)\n",
    "    measureConfig = SingleFrameMeasurementTask.ConfigClass()\n",
    "    measureConfig.plugins.names.add(\"ext_photometryKron_KronFlux\")\n",
    "    measureTask = SingleFrameMeasurementTask(schema=schema, config=measureConfig)\n",
    "    bkgResult = bkgTask.run(exposure)\n",
    "    repairTask.run(exposure)\n",
    "    table = SourceTable.make(schema)\n",
    "    detectResult = detectTask.run(table, exposure)\n",
    "    catalog = detectResult.sources\n",
    "    deblendTask.run(exposure, catalog)\n",
    "    measureTask.run(catalog, exposure)\n",
    "    return catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = processExposure1(exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Add PSF modeling to the end of that function, delegating most of the work to `lsst.pipe.tasks.MeasurePsfTask`.  You may want to use a higher threshold (e.g. 50-sigma) for detection, since PSF modeling should only use bright stars.  Hint: `MeasurePsfTask` requires the catalog you pass it to be contiguous; if it isn't you'll (unfortunately) get a very confusing error message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.pipe.tasks.measurePsf import MeasurePsfTask\n",
    "\n",
    "def processExposure2(exposure):\n",
    "    schema = SourceTable.makeMinimalSchema()\n",
    "    bkgTask = SubtractBackgroundTask()\n",
    "    repairTask = RepairTask()\n",
    "    detectConfig = SourceDetectionTask.ConfigClass()\n",
    "    detectConfig.thresholdValue = 50.0\n",
    "    detectTask = SourceDetectionTask(config=detectConfig, schema=schema)\n",
    "    deblendTask = SourceDeblendTask(schema=schema)\n",
    "    measureConfig = SingleFrameMeasurementTask.ConfigClass()\n",
    "    measureConfig.plugins.names.add(\"ext_photometryKron_KronFlux\")\n",
    "    measureTask = SingleFrameMeasurementTask(schema=schema, config=measureConfig)\n",
    "    psfTask = MeasurePsfTask(schema=schema)\n",
    "    bkgResult = bkgTask.run(exposure)\n",
    "    repairTask.run(exposure)\n",
    "    table = SourceTable.make(schema)\n",
    "    detectResult = detectTask.run(table, exposure)\n",
    "    catalog = detectResult.sources\n",
    "    deblendTask.run(exposure, catalog)\n",
    "    catalog = catalog.copy(deep=True)\n",
    "    measureTask.run(catalog, exposure)\n",
    "    psfTask.run(exposure, catalog)\n",
    "    return catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = processExposure2(exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Make images of the PSF stars, the PSF model at the position of those stars, and the difference between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.image import PARENT\n",
    "\n",
    "stamps = []\n",
    "for record in catalog[catalog[\"calib_psf_used\"]]:\n",
    "    # Get the PSF, convert to a single-precision image (yes, ugly syntax)\n",
    "    psfImage = exposure.getPsf().computeImage(record.getCentroid()).convertF()\n",
    "    # Make a subimage (yes, even uglier syntax).  Last boolean is whether to make a deep copy.\n",
    "    starImage = exposure.image.Factory(exposure.image, psfImage.getBBox(), PARENT, True)\n",
    "    # Divide the star image by the PSF flux (note that this wouldn't work as well if we had applied aperture corrections)\n",
    "    starImage /= record.get(\"slot_PsfFlux_instFlux\")\n",
    "    # Make a deep copy, and subtract to get the residuals\n",
    "    residualImage = starImage.clone()\n",
    "    residualImage -= psfImage\n",
    "    stamps.append((psfImage, starImage, residualImage))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limit = 4\n",
    "nCols = 3\n",
    "nRows = min(len(stamps), limit)\n",
    "index = 1\n",
    "for psfImage, starImage, residualImage in stamps[:nRows]:\n",
    "    vmin = min(psfImage.array.min(), starImage.array.min(), residualImage.array.min())\n",
    "    vmax = max(psfImage.array.max(), starImage.array.max(), residualImage.array.max())\n",
    "    matplotlib.pyplot.subplot(nRows, nCols, index)\n",
    "    display(starImage, vmin=vmin, vmax=vmax)\n",
    "    index += 1\n",
    "    matplotlib.pyplot.subplot(nRows, nCols, index)\n",
    "    display(psfImage, vmin=vmin, vmax=vmax)\n",
    "    index += 1\n",
    "    matplotlib.pyplot.subplot(nRows, nCols, index)\n",
    "    display(residualImage, vmin=vmin, vmax=vmax)\n",
    "    index += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly the PSF model we get without tuning any of the configuration for `MeasurePsfTask` isn't terribly good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5) Add another detect-deblend-measure sequence after PSF modeling at a deeper threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.pipe.tasks.measurePsf import MeasurePsfTask\n",
    "\n",
    "def processExposure3(exposure):\n",
    "    # Construct tasks and schema for the first pass to do PSF estimation\n",
    "    schema1 = SourceTable.makeMinimalSchema()\n",
    "    bkgTask = SubtractBackgroundTask()\n",
    "    repairTask = RepairTask()\n",
    "    detectConfig1 = SourceDetectionTask.ConfigClass()\n",
    "    detectConfig1.thresholdValue = 50.0\n",
    "    detectTask1 = SourceDetectionTask(config=detectConfig1, schema=schema1)\n",
    "    deblendTask1 = SourceDeblendTask(schema=schema1)\n",
    "    measureTask1 = SingleFrameMeasurementTask(schema=schema1)\n",
    "    psfTask = MeasurePsfTask(schema=schema1)\n",
    "    # Construct tasks and schema for the second pass\n",
    "    schema2 = SourceTable.makeMinimalSchema()\n",
    "    detectTask2 = SourceDetectionTask(schema=schema2)\n",
    "    deblendTask2 = SourceDeblendTask(schema=schema2)\n",
    "    measureConfig2 = SingleFrameMeasurementTask.ConfigClass()\n",
    "    measureConfig2.plugins.names.add(\"ext_photometryKron_KronFlux\")\n",
    "    measureTask2 = SingleFrameMeasurementTask(schema=schema2, config=measureConfig2)\n",
    "    # Run the first pass, and do PSF estimation.\n",
    "    bkgResult = bkgTask.run(exposure)\n",
    "    repairTask.run(exposure)\n",
    "    table1 = SourceTable.make(schema1)\n",
    "    detectResult1 = detectTask1.run(table1, exposure)\n",
    "    catalog1 = detectResult1.sources\n",
    "    deblendTask1.run(exposure, catalog1)\n",
    "    catalog1 = catalog1.copy(deep=True)\n",
    "    measureTask1.run(catalog1, exposure)\n",
    "    # Run the second pass\n",
    "    table2 = SourceTable.make(schema2)\n",
    "    detectResult2 = detectTask2.run(table2, exposure)\n",
    "    catalog2 = detectResult2.sources\n",
    "    deblendTask2.run(exposure, catalog2)\n",
    "    measureTask2.run(catalog2, exposure)\n",
    "    return catalog1, catalog2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog1, catalog2 = processExposure3(exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6) Rewrite the function as a class that constructs all of the `Task`s that it use in `__init__` and processes a single `Exposure` with you call its `run` method.  Make sure it will behave properly if `run` is called multiple times wiht different `Exposure` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ProcessExposure(object):\n",
    "    \n",
    "    def __init__(self):\n",
    "        # Construct tasks and schema for the first pass to do PSF estimation\n",
    "        self.schema1 = SourceTable.makeMinimalSchema()\n",
    "        self.bkgTask = SubtractBackgroundTask()\n",
    "        self.repairTask = RepairTask()\n",
    "        detectConfig1 = SourceDetectionTask.ConfigClass()\n",
    "        detectConfig1.thresholdValue = 50.0\n",
    "        self.detectTask1 = SourceDetectionTask(config=detectConfig1, schema=self.schema1)\n",
    "        self.deblendTask1 = SourceDeblendTask(schema=self.schema1)\n",
    "        self.measureTask1 = SingleFrameMeasurementTask(schema=self.schema1)\n",
    "        self.psfTask = MeasurePsfTask(schema=self.schema1)\n",
    "        # Construct tasks and schema for the second pass\n",
    "        self.schema2 = SourceTable.makeMinimalSchema()\n",
    "        self.detectTask2 = SourceDetectionTask(schema=self.schema2)\n",
    "        self.deblendTask2 = SourceDeblendTask(schema=self.schema2)\n",
    "        measureConfig2 = SingleFrameMeasurementTask.ConfigClass()\n",
    "        measureConfig2.plugins.names.add(\"ext_photometryKron_KronFlux\")\n",
    "        self.measureTask2 = SingleFrameMeasurementTask(schema=self.schema2, config=measureConfig2)\n",
    "        \n",
    "    def run(self, exposure):\n",
    "        # Run the first pass, and do PSF estimation.\n",
    "        bkgResult = self.bkgTask.run(exposure)\n",
    "        self.repairTask.run(exposure)\n",
    "        table1 = SourceTable.make(self.schema1)\n",
    "        detectResult1 = self.detectTask1.run(table1, exposure)\n",
    "        catalog1 = detectResult1.sources\n",
    "        self.deblendTask1.run(exposure, catalog1)\n",
    "        catalog1 = catalog1.copy(deep=True)\n",
    "        self.measureTask1.run(catalog1, exposure)\n",
    "        # Run the second pass\n",
    "        table2 = SourceTable.make(self.schema2)\n",
    "        detectResult2 = self.detectTask2.run(table2, exposure)\n",
    "        catalog2 = detectResult2.sources\n",
    "        self.deblendTask2.run(exposure, catalog2)\n",
    "        self.measureTask2.run(catalog2, exposure)\n",
    "        return catalog1, catalog2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure, sub, bbox = regen_exposure()\n",
    "exposure.setPsf(psf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "task = ProcessExposure()\n",
    "catalog1, catalog2 = task.run(exposure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "- "
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}