tutorials/dm_butler_postage_stamps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DC2: Generate Postage Stamps for set of RA, Dec coordinates\n",
    "<br>Owner: **Michael Wood-Vasey** ([@wmwv](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@wmwv))\n",
    "<br>Last Verified to Run: **2018-10-26** (by @yymao)\n",
    "\n",
    "This Notebook demonstrates taking a list of RA, Dec positions and generating \"postage-stamp\" cutout images from the coadded images.  It illustrates using `matplotlib` for simple display and saving of PNGs.  It also introduces using `lsst.afw.display` to use a more DM-product aware visualization to show the detect object footprints and other mask planes.\n",
    "\n",
    "### Learning Objectives:\n",
    "After working through and studying this Notebook you should be able to\n",
    "   1. Generate a postage stamp from a coadd for your chosen RA, Dec in a chosen filter and save as (a) a FITS file and (b) a PNG.\n",
    "   2. Use afw.display to show a postage stamp along with the mask planes from the coadd image, including overlaying the object footprints.\n",
    "   3. [Advanced] Instantiate a Butler object for a repo and read an arbitrary coadd image.  This will require reading through functions carefully and testing individual lines.  \n",
    "   4. [Advanced] Load the PSF model for a coadd and evaluate it at a given RA, Dec.\n",
    "   5. [Expert] Subtract a scaled PSF model from the location of an object on a coadd.\n",
    "   6. [Expert] Look up which tract, patch a given RA, Dec falls in.\n",
    "\n",
    "### Logistics\n",
    "This is intended to be runnable at NERSC through the https://jupyter-dev.nersc.gov interface from a local git clone of https://github.com/LSSTDESC/DC2-analysis in your NERSC directory.  But you can also run it wherever, with appropriate adjustment of the 'repo' location to point to a place where you have a Butler repo will all of the images. Instructions for setting up the proper python kernel can be found here: https://confluence.slac.stanford.edu/x/1_ubDQ\n",
    "\n",
    "Based in part on https://github.com/lsst-com/notebooks/blob/master/postage_stamp.ipynb\n",
    "\n",
    "## Set Up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "\n",
    "from astropy.table import Table\n",
    "\n",
    "import lsst.daf.persistence as dafPersist\n",
    "import lsst.afw.geom as afwGeom\n",
    "import lsst.afw.coord as afwCoord\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.display as afwDisplay\n",
    "\n",
    "from astropy.visualization import ZScaleInterval\n",
    "\n",
    "from desc_dc2_dm_data import REPOS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive backends and NERSC Jupyter vs JupyterLab\n",
    "\n",
    "If you're running on directly in Jupyter and not JupyterLab\n",
    "you can instead use below the interactive figure interface (nbAgg)  \n",
    "`%matplotlib notebook`\n",
    "\n",
    "We don't do that by default in this Notebook so that it can be run under the NERSC JupyterLab environment, which is what DESC specifically called out in the [instructions for using NERSC](https://confluence.slac.stanford.edu/x/1_ubDQ)\n",
    "but I recommend that you use the Jupyter interface to get the best out of this Notebook.\n",
    "\n",
    "If you want to switch from `notebok`<->`inline` you have to restart the Python kernel for this Jupyter Notebook.  I recommend using the \"Kernel->Restart & Clear Output\" menu option above.\n",
    "\n",
    "The distinction between Jupyter and JupyterLab is confusing, so I'll repeat here with the specific examples:\n",
    "1. This is my Jupyter URL and specific URL of this Notebook URL at NERSC:  \n",
    "https://jupyter-dev.nersc.gov/user/wmwv/  \n",
    "https://jupyter-dev.nersc.gov/user/wmwv/notebooks/global/homes/w/wmwv/DC2-analysis/Notebooks/DC2_Postage_Stamps.ipynb\n",
    "\n",
    "2. This is my JupyterLab URL and specific URL of this Notebook at NERSC:  \n",
    "https://jupyter-dev.nersc.gov/user/wmwv/lab  \n",
    "https://jupyter-dev.nersc.gov/user/wmwv/lab/tree/global/homes/w/wmwv/DC2-analysis/Notebooks/DC2_Postage_Stamps.ipynb\n",
    "\n",
    "The interactive 'notebook' (nbAgg) backend will work with #1 but will not work with #2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set plotting to 'inline' figures\n",
    "# If you're running on directly in Jupyter and not JupyterLab\n",
    "# you can instead use the interactive figure interface (nbAgg).\n",
    "%matplotlib inline\n",
    "# %matplotlib notebook\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = (6, 6)\n",
    "zscale = ZScaleInterval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How To Make a Cutout Image\n",
    "\n",
    "The DM butler's `get` method enables postage stamp making, via a bounding box keyword argument. The functions below show how to make this bounding box, and apply it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutout_coadd_ra_dec(butler, ra, dec, filter='r', datasetType='deepCoadd', **kwargs):\n",
    "    \"\"\"\n",
    "    Produce a cutout from coadd from the given butler at the given RA, Dec in decimal degrees.\n",
    "    \n",
    "    Notes\n",
    "    -----\n",
    "    Trivial wrapper around 'cutout_coadd_spherepoint'\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Servant providing access to a data repository\n",
    "    ra: float\n",
    "        Right ascension of the center of the cutout, degrees\n",
    "    dec: float\n",
    "        Declination of the center of the cutout, degrees\n",
    "    filter: string\n",
    "        Filter of the image to load\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \"\"\"\n",
    "    radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)\n",
    "    return cutout_coadd_spherepoint(butler, radec, filter=filter, datasetType=datasetType)\n",
    "    \n",
    "\n",
    "def cutout_coadd_spherepoint(butler, radec, filter='r', datasetType='deepCoadd',\n",
    "                                  skymap=None, cutoutSideLength=51, **kwargs):\n",
    "    \"\"\"\n",
    "    Produce a cutout from a coadd at the given afw SpherePoint radec position.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Servant providing access to a data repository\n",
    "    radec: lsst.afw.geom.SpherePoint \n",
    "        Coordinates of the center of the cutout.\n",
    "    filter: string \n",
    "        Filter of the image to load\n",
    "    datasetType: string ['deepCoadd']  \n",
    "        Which type of coadd to load.  Doesn't support 'calexp'\n",
    "    skymap: lsst.afw.skyMap.SkyMap [optional] \n",
    "        Pass in to avoid the Butler read.  Useful if you have lots of them.\n",
    "    cutoutSideLength: float [optional] \n",
    "        Side of the cutout region in pixels.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \"\"\"\n",
    "    cutoutSize = afwGeom.ExtentI(cutoutSideLength, cutoutSideLength)\n",
    "\n",
    "    if skymap is None:\n",
    "        skymap = butler.get(\"%s_skyMap\" % datasetType)\n",
    "    \n",
    "    # Look up the tract, patch for the RA, Dec\n",
    "    tractInfo = skymap.findTract(radec)\n",
    "    patchInfo = tractInfo.findPatch(radec)\n",
    "    xy = afwGeom.PointI(tractInfo.getWcs().skyToPixel(radec))\n",
    "    bbox = afwGeom.BoxI(xy - cutoutSize//2, cutoutSize)\n",
    "\n",
    "    coaddId = {'tract': tractInfo.getId(), 'patch': \"%d,%d\" % patchInfo.getIndex(), 'filter': filter}\n",
    "    \n",
    "    cutout_image = butler.get(datasetType+'_sub', bbox=bbox, immediate=True, dataId=coaddId)\n",
    "    \n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_cutout_image(butler, ra, dec, filter='r', vmin=None, vmax=None, label=None,\n",
    "                      show=True, saveplot=False, savefits=False,\n",
    "                      datasetType='deepCoadd'):\n",
    "    \"\"\"\n",
    "    Generate and optionally display and save a postage stamp for a given RA, Dec.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Servant providing access to a data repository\n",
    "    ra: float\n",
    "        Right ascension of the center of the cutout, degrees\n",
    "    dec: float\n",
    "        Declination of the center of the cutout, degrees\n",
    "    filter: string \n",
    "        Filter of the image to load\n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Uses matplotlib to generate stamps.  Saves FITS file if requested.\n",
    "    \"\"\"\n",
    "\n",
    "    cutout_image = cutout_coadd_ra_dec(butler, ra, dec, filter=filter, datasetType='deepCoadd')\n",
    "    if savefits:\n",
    "        if isinstance(savefits, str):\n",
    "            filename = savefits\n",
    "        else:\n",
    "            filename = 'postage-stamp.fits'\n",
    "        cutout_image.writeFits(filename)\n",
    "    \n",
    "    radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)\n",
    "    xy = cutout_image.getWcs().skyToPixel(radec)\n",
    "    \n",
    "    if vmin is None or vmax is None:\n",
    "        vmin, vmax = zscale.get_limits(cutout_image.image.array)\n",
    "\n",
    "    plt.imshow(cutout_image.image.array, vmin=vmin, vmax=vmax, cmap='binary_r', origin='lower')\n",
    "    plt.colorbar()\n",
    "    plt.scatter(xy.getX() - cutout_image.getX0(), xy.getY() - cutout_image.getY0(),\n",
    "                color='none', edgecolor='red', marker='o', s=200)\n",
    "    if label is not None:\n",
    "        plt.title(label)\n",
    "    if saveplot:\n",
    "        if isinstance(saveplot, str):\n",
    "            filename = saveplot\n",
    "        else:\n",
    "            filename = 'postage-stamp.png'\n",
    "        plt.savefig(filename)\n",
    "    if show:\n",
    "        plt.show()\n",
    "\n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making Some Postage Stamps\n",
    "\n",
    "We'll use a list of 20 targets, pre-selected and checked into the Notebooks folder, and have the butler make postage stamps from the Run 1.1 coadd images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run = '1.1p'\n",
    "repo = REPOS[run]\n",
    "butler = dafPersist.Butler(repo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter = 'r'\n",
    "coord_file = 'assets/id_ra_dec_mid_mag_{}_{}.txt'.format(run, filter)\n",
    "id_ra_dec = Table.read(coord_file, format='ascii')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot just one\n",
    "i = 3\n",
    "first = id_ra_dec[i]\n",
    "ra, dec = first['RA'], first['DEC']\n",
    "frame = 0\n",
    "plt.figure(frame)\n",
    "cutout = make_cutout_image(butler, ra, dec, filter=filter, label=\"Object ID: %d\" % id_ra_dec[i]['ID'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use `notebook` as the backend, we can interact with the figure and read off pixel coordinates and counts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can loop over this to create postage-stamps for each entry in our catalog.  To keep the Notebooks folder tidy, we'll make and use a subfolder.\n",
    "\n",
    "We use `show=False` because we don't necessarily want to display each of the 20 figures.  But even with `show=False` we will get a figure displayed in interactive `%matplotlib notebook` mode because under that backend, one is always displaying the figure without waiting for a `plt.show` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subfolder = 'Stamps'\n",
    "if not os.path.exists(subfolder):\n",
    "    os.mkdir(subfolder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = 1\n",
    "show = False\n",
    "datasetType = 'deepCoadd'\n",
    "vmin, vmax = -0.75, +1.5  # Fix the vmin, vmax to make it easier to compare across postage stamps.\n",
    "\n",
    "for objectId, ra, dec in id_ra_dec:\n",
    "    plt.figure(frame)\n",
    "    plt.clf()\n",
    "    basename = \"Stamps/%s_%s_%s\" % (datasetType, objectId, filter)\n",
    "    saveplot = \"%s.png\" % basename\n",
    "    savefits = \"%s.fits\" % basename\n",
    "    make_cutout_image(butler, ra, dec, filter=filter, vmin=vmin, vmax=vmax,\n",
    "                      datasetType=datasetType,\n",
    "                      label=\"Object ID: %d\" % objectId,\n",
    "                      show=show, saveplot=saveplot, savefits=savefits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## AFW Display\n",
    "\n",
    "But there's additional information available in the Exposure object that we can access using LSST DM tools that are aware of these.  Specifically, we'll use `lsst.afw.display` to expose the mask planes in a cutout image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def display_cutout_image(butler, ra, dec, vmin=None, vmax=None, label=None,\n",
    "                      frame=None, display=None, backend='matplotlib',\n",
    "                      show=True, saveplot=False, savefits=False,\n",
    "                      old_matplotlib = False,\n",
    "                      datasetType='deepCoadd'):\n",
    "    \"\"\"\n",
    "    Display a postage stamp for a given RA, Dec using LSST lsst.afw.display.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    backend: string\n",
    "        Backend can be anything that lsst.afw.display and your configuration supports: \n",
    "        e.g. matplotlib, ds9, ginga, firefly.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \n",
    "    Notes\n",
    "    -----\n",
    "    Parameters are the same as for make_cutout_image, except for the backend.\n",
    "    You definitely have the matplotlib backend.\n",
    "    ds9, ginga, and firefly can be set up but are non-trivial on the scale of a simple Notebook.\n",
    "    \"\"\"\n",
    "    cutout_image = cutout_coadd_ra_dec(butler, ra, dec, filter='r', datasetType='deepCoadd')\n",
    "    if savefits:\n",
    "        if isinstance(savefits, str):\n",
    "            filename = savefits\n",
    "        else:\n",
    "            filename = 'postage-stamp.fits'\n",
    "        cutout_image.writeFits(filename)\n",
    "    \n",
    "    if display is None:\n",
    "        display = afwDisplay.Display(frame=frame, backend=backend)\n",
    "\n",
    "    radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)\n",
    "    xy = cutout_image.getWcs().skyToPixel(radec)\n",
    "    \n",
    "    display.mtv(cutout_image)\n",
    "    display.scale(\"asinh\", \"zscale\")\n",
    "    display.dot('o', xy.getX(), xy.getY(), ctype='red')\n",
    "    display.show_colorbar()\n",
    "    \n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display just one\n",
    "i = 3\n",
    "first = id_ra_dec[i]\n",
    "ra, dec = first['RA'], first['DEC']\n",
    "\n",
    "# We start specifying figure numbers here and below to explicitly control which figure we're plotting to\n",
    "# If we use the interactive `%matplotlib notebook` (nbAgg),\n",
    "# then figures stay active even when you exit the cell.\n",
    "# It can then be confusing to figure out which figure is active.\n",
    "# (matplotlib.pyplot doesn't call these 'frame's, but afwDisplay does)\n",
    "frame = 2\n",
    "plt.figure(frame)\n",
    "cutout = display_cutout_image(butler, ra, dec, frame=frame,\n",
    "                              label=\"Object ID: %d\" % id_ra_dec[i]['ID'])\n",
    "\n",
    "# After we have the displayed image we can interact with it a bit through matplotlib\n",
    "plt.savefig('afwDisplay.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here when we browse around we get both x,y coordinates (RA, Dec) = $(\\alpha, \\delta)$ the pixel counts as well as the list of named mask bits that are set for the pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using mask planes\n",
    "\n",
    "Along with the pixels, the coadd image also has a set of associated \"mask\" planes. \"Mask\" means has an identified property -- it doesn't necessarily mean \"bad\".  In particular, the \"DETECTED\" mask plane means the measurement identified an object here.  \"Plane\" here refers to a specific bit in the mask.\n",
    "\n",
    "Specifically, the blue overlays above are called the \"footprints\" for the observations.  These are the pixels identified by the pipeline as \"belonging\" to the object: pixels where the object contributes detectable number of photons above the sky background.\n",
    "\n",
    "The colors are configurable, but the default values are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = 1\n",
    "plt.figure(frame)\n",
    "display = afwDisplay.Display(frame=frame, backend='matplotlib')\n",
    "for maskName, maskBit in cutout.mask.getMaskPlaneDict().items():\n",
    "    print('{}: {}'.format(maskName, display.getMaskPlaneColor(maskName)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing `plt.imshow` and `afwDisplay`:\n",
    "* X,Y origin  \n",
    "    Note that `afwDisplay` displays images with the lower-left as the origin, which is the x,y convention that most  astronomers are using to thinking in.  This is opposite to the default vertical orientation in `matplotlib.pyplot.imshow`, thus we passed `origin='lower'` in our call above to `imshow`. `afwDisplay` also provides the original pixel coordinates of the tract+patch coadd image, which is potentially useful.\n",
    "* Colormap  \n",
    "    This is an arbitrary choice, but `afwDisplay` displays pixels with large values in lighter shades. We match this with the `binary_r` colormap in `matplotlib`.\n",
    "* Colorbar \n",
    "  We have added the colorbar explicitly with a call to `plt.colorbar`.\n",
    "* Saving output\n",
    "  Again, we've used the `matplotlib` interface to save the resulting figure to a PNG.  This is not presently possible directly from the afwDisplay object.\n",
    "\n",
    "Take a look at https://pipelines.lsst.io/v/DM-11077/getting-started/display.html for a few more details on afwDisplay. You can also look at [the source code](https://github.com/lsst/display_matplotlib/blob/master/python/lsst/display/matplotlib/matplotlib.py) for the matplotlib backend to afwDisplay."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the PSF Model\n",
    "\n",
    "The cutout image object retains full information from the image, including the PSF model.\n",
    "\n",
    "The PSF model is accessible as a function object that can be evaluated a specific locations.\n",
    "\n",
    "In addition, the coadd catalog saves the xx, xy, yy moments of the PSF model at the location of photometered objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We still have our cutout image object from above\n",
    "psf = cutout.getPsf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(psf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The repetition in the object name is a consequence of the pybind11 mapping.  It's really a `lsst.meas.algorithms.coaddPsf.CoaddPsf` object.  The C++ documentation unfortunately doesn't come through nicely through the wrapping to Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(psf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look up the documentation in the LSST DM Stack Doxygen pages.  If we assume (correctly for these purposes) that a CoaddPsf is just like a regular PSF, then we can look here:\n",
    "http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1afw_1_1detection_1_1_psf.html\n",
    "\n",
    "where we will see that there are several relevant functions:\n",
    "`computeImage`: Return an Image of the PSF, in a form that can be compared directly with star images.\n",
    "`computeKernelImage`:  Return an Image of the PSF, in a form suitable for convolution.\n",
    "`computeShape`: Compute the ellipse corresponding to the second moments of the Psf.\n",
    "\n",
    "and then some interesting things you might not have expected to even be available:\n",
    "`getAverageColor`: Return the average Color of the stars used to construct the Psf.\n",
    "`getAveragePosition`: Return the average position of the stars used to construct the Psf. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the documentation is hard to read from a Python perspective.  What kind of arguments does the following function actually want and how do I create those objects?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "computeImage(...) from builtins.PyCapsule\n",
    "        computeImage(self: lsst.afw.detection._psf.Psf, position: lsst.afw.geom.coordinates.coordinates.Point2D=Point2D(nan, nan), color: lsst.afw.image.color.Color=<lsst.afw.image.color.Color object at 0x2b8f57b38928>, owner: lsst.afw.detection._psf.ImageOwnerEnum=ImageOwnerEnum.COPY) -> lsst.afw.image.image.image.ImageD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examples are easiest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's use the same first row to get an RA, Dec and translate that to an x, y on the image plane.\n",
    "first = id_ra_dec[i]\n",
    "ra, dec = first['RA'], first['DEC']\n",
    "radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)\n",
    "xy = cutout.getWcs().skyToPixel(radec)  # returns a Point2D\n",
    "\n",
    "kernel_image = psf.computeImage(xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object we get back is an `lsst.afw.image.image.ImageD` -- An `Image` in double precision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's display the kernel image we made:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(kernel_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, the postage stamp is a full `lsst.afw.image.exposure.ExposureF`.  An `Exposure` contains a `MaskedImage` (an image, a mask, and variance plane), along with an optional SkyWCS.\n",
    "http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1afw_1_1image_1_1_exposure.html#details\n",
    "\n",
    "An `Exposure` can also generate `subExposures`, which are `Exposure`s of a subset of the pixels.  Our postage stamp cutout is thus an `Exposure`.  The additional information is what allowed `afwDisplay` to show the image with the original pixel coordinates.\n",
    "\n",
    "If you read a `raw`, `calexp`, or `deepCoadd` or similar image product through the Butler, you will get an `Exposure`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(cutout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = 2\n",
    "plt.figure(frame)\n",
    "display = afwDisplay.Display(frame=frame, backend='matplotlib')\n",
    "display.scale(\"linear\", \"minmax\")\n",
    "display.mtv(kernel_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we see an image of the PSF centered at the x, y position corresponding to the RA, Dec position of the object.  This is a convenient representation to use to subtract.\n",
    "\n",
    "Note that while the cutout was 51x51 pixels (the default `cutoutSideLength` in `cutout_coadd_spherepoint` above), the PSF is 61x61 pixels.  That's the default size of specifying the full PSF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Cutout shape: \", np.shape(cutout.image.array))\n",
    "print(\"PSF shape: \", np.shape(kernel_image.array))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sum(kernel_image.array))  # Kernel is normalized to total sum of 1 over the size of the PSF."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next scale our PSF model to the detected object in the coadd and illustrate using a 1D slice.  We have to account for the different sizes of the PSF model and the cutout, including the finite fraction of PSF flux included in the smaller cutout  The scaling is the the total sum along the slice in this overly simple model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = 3\n",
    "plt.figure(frame)\n",
    "cutout_nx, cutout_ny = np.shape(cutout.image.array)\n",
    "kernel_nx, kernel_ny = np.shape(kernel_image.array)\n",
    "\n",
    "cutout_slice = cutout.image.array[cutout_nx//2, :]\n",
    "kernel_slice = kernel_image.array[kernel_nx//2, :]\n",
    "\n",
    "# Normalize PSF sum to cutout sum in this slice\n",
    "scaling = np.sum(cutout_slice)/np.sum(kernel_slice)\n",
    "\n",
    "plt.step(np.arange(cutout_nx) - cutout_nx/2, cutout_slice,\n",
    "         label='Postage Stamp slice',\n",
    "         linewidth=2)\n",
    "plt.step(np.arange(kernel_nx) - kernel_nx/2, scaling * kernel_slice,\n",
    "         label='PSF scaled by %.2f' % scaling,\n",
    "         linewidth=2)\n",
    "\n",
    "# Assume that the PSF kernel size is *larger* than the cutout size.\n",
    "offset_nx = kernel_nx - cutout_nx\n",
    "offset_ny = kernel_ny - cutout_ny\n",
    "\n",
    "plt.step(np.arange(cutout_nx) - cutout_nx/2,\n",
    "         cutout_slice - (scaling * kernel_slice[offset_nx//2:-offset_nx//2]),\n",
    "         label='Postage stamp - scaled PSF',\n",
    "         linestyle='--', color='black')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('counts in pixel')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subtract a scaled version of the PSF from the image array.\n",
    "Note that we have to select out just the section of the PSF that's within the cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutout_minus_psf = cutout.image.array - \\\n",
    "    scaling * kernel_image.array[offset_nx//2:-offset_nx//2,\n",
    "                                 offset_ny//2:-offset_ny//2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create an `Exposure` with the pixel values of subtracting the PSF model we create a copy of the postage stamp and set the image data array to the residual value we just calculated above.\n",
    "We keep the WCS, the pixel mapping of this subExposure, and the mask and variance plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residual = cutout.clone()\n",
    "residual.setImage(afwImage.ImageF(afwImage.ImageF(cutout_minus_psf.astype(np.float32))))\n",
    "\n",
    "frame = 4\n",
    "plt.figure(frame)\n",
    "display = afwDisplay.Display(frame=frame,backend='matplotlib')\n",
    "display.scale(\"linear\", \"minmax\")\n",
    "display.mtv(cutout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = 5\n",
    "fig = plt.figure(frame)\n",
    "\n",
    "display = afwDisplay.Display(frame=frame, backend='matplotlib')\n",
    "\n",
    "display.scale(\"linear\", \"minmax\")\n",
    "display.mtv(residual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the source is gone.  The actual fitting for the PSF flux is a bit more complicated, but the above demonstrates the basic process of evaluating a PSF at a location in the image and creating a new `Image` or `Exposure` object with a PSF model subtracted out."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-stack",
   "language": "python",
   "name": "desc-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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}