workshops/lyon2018/intro-process-ccd.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.  Some of the things we'll be doing are more commonly done on the command-line, via executable scripts the stack also provides.  More complete tutorials for data processing and command-line functionality can be found in the [Science Pipelines documentation](https://pipelines.lsst.io/getting-started/index.html) and [DM Tech Note 23](https://dmtn-023.lsst.io/).\n",
    "\n",
    "More notebook examples can be found here: https://github.com/RobertLuptonTheGood/notebooks/tree/master/Demos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Repository Setup\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*.  The butler operates on *data repositories*.  For this tutorial we have setup a data repository on lsp-demo, at `/project/ci_hsc_small` read-only space. This data repository contains some raw data, master calibration files, transmission curves data, and an external reference catalog.  All of these are from a self-contained test dataset we call [ci_hsc](https://github.com/lsst/ci_hsc).  The full ci_hsc dataset includes just enough data to run the full (current) LSST pipeline, which extends through processing coadds from multiple bands together.  In this tutorial we'll focus on processing an individual image, and that's all this particular subset will support.  We also won't go into the details of how to build master calibration files or reference catalogs here.\n",
    "\n",
    "First, copy the data repository to your home space, where you have full permission to create and manage data.  Open a terminal window and do:\n",
    "\n",
    "    cp -r /project/ci_hsc_small/* $HOME/WORK\n",
    "\n",
    "In this tutorial `$HOME/WORK` is your working data repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "WORK_DIR = os.path.join(os.environ['HOME'], \"WORK\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instrument Signature Removal and Command-Line Tasks\n",
    "\n",
    "Before we can start doing *interesting* things, we need some minimally processed images (i.e. flat-fielded, bias-corrected, etc).  Because the HSC team has spent a lot of time characterizing the instrument, we really want to run this step with the default configuration they've provided.  That's also much actually easier to do from the command-line, and while we *could* do it from Python, that'd involve a lot of little irrelevant workarounds we'd rather not get bogged down in.\n",
    "\n",
    "ISR is implemented as a subclass of `lsst.pipe.base.Task`.  Nearly all of our high-level algorithms are implemented as `Task`s, which are essentially just callable objects that can be composed (a high-level `Task` can hold one or more lower-level \"subtasks\", to which it can delegate work) and configured (every task takes an instance of a configuration class that controls what it does in detail).  ISR is actually a `CmdLineTask`, a special kind of task that can be run from the command-line and use the data butler for all of its inputs and outputs (regular `Task`s generally do not use the butler directly).  Unlike virtually every other algorithm, there is a different ISR `Task` for each major camera (though there's also a simple default one), reflecting the specialized processing that's needed at this level.\n",
    "\n",
    "But (for uninteresting, historical reasons), it's not currently possible to run `IsrTask` from the command-line.  Instead, what we can do is run a parent `CmdLineTask`, `lsst.pipe.tasks.ProcessCcdTask`, which will run `IsrTask` as well as a few other steps.  By default, it doesn't actually save the image directly after ISR is run - it performs a few more operations first, and then saves that image.  But we can tell it to do so by modifying the tasks' configuration when we.\n",
    "\n",
    "In a terminal session, load the LSST environment and setup the software stack by doing:\n",
    "\n",
    "    source /opt/lsst/software/stack/loadLSST.bash\n",
    "    setup lsst_distrib\n",
    "\n",
    "The full command-line for running `ProcessCcdTask` is below:\n",
    "\n",
    "    processCcd.py $HOME/WORK --rerun isr --id visit=903334 ccd=16 --config isr.doWrite=True\n",
    "\n",
    "Alternatively, we can use a Jupyter trick (the `%%script bash` line) to do that within the notebook: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%script bash\n",
    "processCcd.py $HOME/WORK --rerun isr --id visit=903334 ccd=16 --config isr.doWrite=True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few features of this command-line that bear explaining:\n",
    " - We run `processCcd.py`, not `ProcessCcdTask`.  There's a similar driver script for all `CmdLineTask`s, with the name formed by making the first word lowercase and removing the `Task` suffix.  These are added to your `PATH` when you set up the LSST package in which they're defined (which happens automatically in the JupyterLab environment).\n",
    " - The first argument to any command-line task is the path an input data repository.\n",
    " - We've used the `--rerun` argument to set the location of the output repository, in this case `$HOME/WORK/rerun/isr`.  You can also use `--output` to set the path more directly, but we recommend `--rerun` because it enforces a nice convention for where to put outputs that helps with discoverability.\n",
    " - The `--id` argument sets the *data ID(s)* to be processed, in this case a single CCD from a single visit.  All `CmdLineTasks` share a fairly sophisticated syntax for expressions that match multiple data IDs, which you can learn more about by running any `CmdLineTask` with `--help`.\n",
    " - We've overridden a configuration value with the `--config` option, in this case to make sure the just-after-ISR image file is written.  Running a `CmdLineTask` automatically also includes applying configuation overrides that customize the task for the kind of data you're processing (i.e. which camera it comes from), and that's how the task knows to run the custom ISR task for HSC, rather than the generic default.  You can see all of the config options for a `CmdLineTask` by running with `--show config`, though the results can be a bit overwhelming.\n",
    " \n",
    "The rest of this tutorial is focused on using LSST software as a Python library, so this will be the last thing we run from the command-line.  Again, for more information about how to run LSST's existing processing scripts from the command-line, check out [DM Tech Note 23](https://dmtn-023.lsst.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Access with Butler\n",
    "\n",
    "The outputs of `CmdLineTasks`, like their inputs, are organized into data repositories, which are managed by an object called `Butler`.  To retrieve a dataset from the `Butler`, we start by constructing one pointing to the output repository from the processing run (which is now an input repository for this `Butler`, which won't have an output repository since we won't be writing any more files):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.daf.persistence import Butler\n",
    "butler = Butler(inputs=os.path.join(WORK_DIR, \"rerun/isr\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then call `get` with the name and data ID of the dataset.  The name of the image that's saved directly after ISR is `postISRCCD`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure = butler.get(\"postISRCCD\", visit=903334, ccd=16)"
   ]
  },
  {
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.afw.geom import Box2D, Box2I, Point2I, Extent2I\n",
    "from lsst.afw.image import Exposure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell (and the one below) to re-load the post-ISR Exposure from disk after\n",
    "# modifying it.\n",
    "exposure = butler.get(\"postISRCCD\", visit=903334, ccd=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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=175, vmax=300, 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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SubtractBackgroundTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bkgConfig = SubtractBackgroundTask.ConfigClass()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "help(bkgConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SubtractBackgroundTask.ConfigClass.algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bkgTask = SubtractBackgroundTask(config=bkgConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": [
    "## 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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SingleGaussianPsf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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(x, y)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.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(x, y)`.  That will shift the image returned by `computeKernelImage(x, y)` 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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.pipe.tasks.repair import RepairTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "repairTask = RepairTask()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.meas.algorithms import SourceDetectionTask\n",
    "from lsst.afw.table import SourceTable, SourceCatalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "schema = SourceTable.makeMinimalSchema()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "detectTask = SourceDetectionTask(schema=schema)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.meas.deblender import SourceDeblendTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deblendTask = SourceDeblendTask(schema=schema)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "catalog = detectResult.sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lsst.meas.base import SingleFrameMeasurementTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Configure the new measurement to run\n",
    "measureConfig.plugins.names.add(\"ext_photometryKron_KronFlux\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "measureTask = SingleFrameMeasurementTask(schema=schema, config=measureConfig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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",
    "    patch = matplotlib.patches.Ellipse((record.getX(), record.getY()),\n",
    "                                       axes.getA(), axes.getB(), axes.getTheta() * 180.0 / np.pi,\n",
    "                                      fill=False, edgecolor=\"green\")\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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "catalog = catalog.copy(deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psfFlux = catalog[\"base_PsfFlux_flux\"]"
   ]
  },
  {
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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",
    "0. Make a scatter plot of Kron Flux vs. Psf Flux.\n",
    "\n",
    "1. 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.\n",
    "\n",
    "2. 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.\n",
    "\n",
    "3. Make images of the PSF stars, the PSF model at the position of those stars, and the difference between them.\n",
    "\n",
    "4. Add another detect-deblend-measure sequence after PSF modeling at a deeper threshold.\n",
    "\n",
    "5. 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": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "- "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST_Stack (Python 3)",
   "language": "python",
   "name": "lsst_stack"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}