Basic DM AFW Introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Introduction to the DM AFW packages\n",
    "#### Written by Chris Walter of Duke University for the LSST Dark Energy Science Collaboration. Last Updated 5/2015 by Will Dawson using LSST Stack v10_0.\n",
    "\n",
    "This IPython notebook is an introduction to the LSST DM stack with a focus on the AFW (applications framework) package used for dealing with image data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Don't forget to the load LSST setup files and then \"setup\" before you start.\n",
    "For this notebook running 'setup pipe_tasks' should load what you need. So before runing `ipython notebook` run the following in the terminal:\n",
    "\n",
    "```\n",
    "cd \\path\\to\\lsst\\\n",
    "source loadLSST.bash\n",
    "setup pipe_tasks -T v10_0\n",
    "ipython notebook\n",
    "```\n",
    "\n",
    "If you followed the default LSST DM pipeline installation instructions then \\path\\to\\lsst\\ is ~\\lsst. Change the suffix of the `source loadLSST.bash` command to match your shell.\n",
    "\n",
    "Start by importing the packages we are going to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import lsst.afw.math        as math\n",
    "import lsst.afw.table       as afwTable\n",
    "import lsst.afw.image       as afwImg\n",
    "import lsst.afw.detection   as afwDetect\n",
    "import lsst.meas.algorithms as measAlg\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib.colors  import LogNorm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exposure Creation and Exploration\n",
    "Now create an exposure by reading a PhoSim output file (you should have created this seperately and copied the eimage file to somewhere you can use it).  A exposure contains a MaskedImage (image, mask and variance image), Meta-data and optional WCS and PSF information. The mask and variances will be empty. The Psf returns nothing now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "exposure    = afwImg.ExposureF('lsst_e_99999999_f2_R22_S11_E000.fits.gz')\n",
    "\n",
    "maskedImage = exposure.getMaskedImage()\n",
    "\n",
    "# These three are held in the maskedImage\n",
    "image       = maskedImage.getImage()\n",
    "mask        = maskedImage.getMask()\n",
    "variance    = maskedImage.getVariance()\n",
    "\n",
    "wcs         = exposure.getWcs()\n",
    "metaData    = exposure.getMetadata()\n",
    "psf         = exposure.getPsf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the CCD data (there should be a star in the middle) and also print out the mask and variance data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAEACAYAAAB/KfmzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADpJJREFUeJzt3G+s3mV9x/H3hz+NOJmkMWn50wUelIQuWzBd6LLhKJtj\nxSxAlgX0gSHa7AluGE2c1Cyjj9SYLI49gD2YjsK2mmZmDCMqlVA1S6T7Q6VSK5TZxJ7Zw4zMPzHL\nSvzuwX0d++NwTs9pr/u+T0/7fiUnvX7X/ftd39/F6f3p7y+pKiSpxwUrvQOSVj+DRFI3g0RSN4NE\nUjeDRFI3g0RSt6kHSZJtSQ4neTHJh6ddX9L4ZZrPkSS5EPg28HZgBvhX4F1V9a2p7YSksZv2EckN\nwJGqOlpVJ4DPALdPeR8kjdm0g+RK4LuD5WOtT9IqNu0g8Xl86Rx00ZTrzQAbBssbGB2V/FwSw0aa\nkqrKuAaa2g+j4HoJuBpYAxwArpu3Tk1xf3aei7Wcm7WWWa/GNdZUj0iq6tUkfwx8CbgQ+FR5x0Za\n9aZ9akNVfQH4wrTrSpqc8/3J1n3naK1p15tmrWnXO1drjdVUH0hbjiRV47oAJGlR4/yune9HJJLG\nwCCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTN\nIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0g\nkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR1M0gkdTNIJHUzSCR\n1G3JIEny6SSzSQ4O+tYm2ZvkhSRPJrls8NmOJC8mOZzklkH/5iQH22cPjH8qklbKco5I/hbYNq/v\nPmBvVV0LPNWWSbIJuAvY1LZ5MEnaNg8B26tqI7AxyfwxJa1SSwZJVX0NeGVe923ArtbeBdzR2rcD\nu6vqRFUdBY4AW5JcDlxaVfvbeo8MtpG0yp3pNZJ1VTXb2rPAuta+Ajg2WO8YcOUC/TOtX9I54KLe\nAaqqktQ4dmZOkp2DxX1VtW+c40vnoyRbga2TGPtMg2Q2yfqqOt5OW15u/TPAhsF6VzE6Eplp7WH/\nzGKDV9XOM9wvSYto/yDvm1tOcv+4xj7TU5vHgbtb+27gsUH/O5OsSXINsBHYX1XHgR8l2dIuvr57\nsI2kVW7JI5Iku4GbgLck+S7w58DHgT1JtgNHgTsBqupQkj3AIeBV4J6qmjvtuQd4GLgEeKKqvjje\nqUhaKTn5PT87JKmqytJrSuoxzu+aT7ZK6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQSOpmkEjq\nZpBI6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQSOpm\nkEjqZpBI6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQSOpmkEjqZpBI6maQ\nSOpmkEjqZpBI6maQSOpmkGjqkmu2rvQ+aLxSVSu9D6+RpKoqK70f0rlunN81j0gkdTNIJHUzSCR1\nWzJIkmxI8nSS55N8M8m9rX9tkr1JXkjyZJLLBtvsSPJiksNJbhn0b05ysH32wGSmJGnalnNEcgL4\nQFX9MvDrwPuSXAfcB+ytqmuBp9oySTYBdwGbgG3Ag0nmLug8BGyvqo3AxiTbxjobSStiySCpquNV\ndaC1fwJ8C7gSuA3Y1VbbBdzR2rcDu6vqRFUdBY4AW5JcDlxaVfvbeo8MtpG0ip3WNZIkVwNvBZ4B\n1lXVbPtoFljX2lcAxwabHWMUPPP7Z1q/pFXuouWumORNwGeB91fVj0+erUBVVZKxPZCSZOdgcV9V\n7RvX2NL5KslWYOskxl5WkCS5mFGIPFpVj7Xu2STrq+p4O215ufXPABsGm1/F6EhkprWH/TML1auq\nncuegVZc1vxz8eZLv8r3//A9VT/4z5XeHy2s/YO8b245yf3jGns5d20CfAo4VFV/OfjoceDu1r4b\neGzQ/84ka5JcA2wE9lfVceBHSba0Md892Ear2Yk/exfff+iv4bo/WOld0cpY8hH5JDcCXwWeA+ZW\n3gHsB/YAvwQcBe6sqv9p23wEeC/wKqNToS+1/s3Aw8AlwBNVde8C9XxEXpqCcX7XfNdGOk/5ro2k\ns4pBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKp\nm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmb\nQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKpm0EiqZtB\nIqmbQSKpm0EiqdspgyTJG5I8k+RAkkNJPtb61ybZm+SFJE8muWywzY4kLyY5nOSWQf/mJAfbZw9M\nbkqSpu2UQVJV/wvcXFXXA78K3JzkRuA+YG9VXQs81ZZJsgm4C9gEbAMeTJI23EPA9qraCGxMsm0S\nE5I0fUue2lTVT1tzDXAh8ApwG7Cr9e8C7mjt24HdVXWiqo4CR4AtSS4HLq2q/W29RwbbSFrllgyS\nJBckOQDMAk9X1fPAuqqabavMAuta+wrg2GDzY8CVC/TPtH5J54CLllqhqn4GXJ/kzcCXktw87/NK\nUuPcqSQ7B4v7qmrfOMeXzkdJtgJbJzH2kkEyp6p+mOTzwGZgNsn6qjreTltebqvNABsGm13F6Ehk\nprWH/TOnqLVzufslaXnaP8j75paT3D+usZe6a/OWuTsySS4Bfhd4FngcuLutdjfwWGs/DrwzyZok\n1wAbgf1VdRz4UZIt7eLruwfbSFrlljoiuRzYleQCRqHzaFU9leRZYE+S7cBR4E6AqjqUZA9wCHgV\nuKeq5k577gEeBi4BnqiqL457MpJWRk5+z88OSaqqsvSaknqM87vmk62SuhkkkroZJJK6GSSSuhkk\nkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSS\nuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6\nGSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6GSSSuhkkkroZJJK6LStIklyY5Nkk\nn2vLa5PsTfJCkieTXDZYd0eSF5McTnLLoH9zkoPtswfGPxVJK2W5RyTvBw4B1ZbvA/ZW1bXAU22Z\nJJuAu4BNwDbgwSRp2zwEbK+qjcDGJNvGMwVJK23JIElyFfAO4G+AuVC4DdjV2ruAO1r7dmB3VZ2o\nqqPAEWBLksuBS6tqf1vvkcE2kla55RyRfBL4EPCzQd+6qppt7VlgXWtfARwbrHcMuHKB/pnWL+kc\ncMogSfL7wMtV9Swnj0Zeo6qKk6c8ks5DFy3x+W8AtyV5B/AG4BeTPArMJllfVcfbacvLbf0ZYMNg\n+6sYHYnMtPawf2axokl2Dhb3VdW+ZczltCXZOqmxV7LWtOs5t9VRK8lWYOtEBq+qZf0ANwGfa+1P\nAB9u7fuAj7f2JuAAsAa4BngJSPvsGWALoyObJ4Bti9Sp5e5T7w+w81ys5dystcx6Na6xljoieV3u\ntD8/DuxJsh04CtzZ9upQkj2M7vC8CtxTbY+Be4CHgUuAJ6rqi6dZW9JZatlBUlVfAb7S2j8A3r7I\neh8FPrpA/78Dv3JmuynpbJaTBwxnhyRn1w5J57CqWvAmyuk664JE0urjuzaSuhkkkvpN83bTErei\ntgGHgRdpt5bHMOZR4DngWWB/61sL7AVeAJ4ELhusv6PVPwzcsozxP83oyd6Dg77THh/YDBxsnz1w\nGrV2MnpO59n2c+uYam0AngaeB74J3DvhuS1Wb+zzY/Q81DOMHlM4BHxsUnM7Ra2J/N4G617Yxv3c\npOb2uprTCIllfCEvZPReztXAxe0//HVjGPc7wNp5fZ8A/rS1P8zrn4G5uO3HEeCCJcZ/G/BWXvvl\nPp3x565R7QduaO0Fn7FZpNb9wAcXWLe31nrg+tZ+E/Bt4LoJzm2xepOa3xvbnxcBXwdunODcFqo1\nkXkNxvkg8PfA45P8Ozn8OVtObW4AjlTV0ao6AXyG0QuA4zD/qvTpvHB4w6kGrqqvAa90jL/sFxoX\nqbXQ/MZR63hVHWjtnwDfYvRu1KTmtli9Sc3vp625htE/Yq9McG4L1ZrIvGDlXrI9W4LkSuC7g+W5\nl/16FfDlJP+W5I9a3+m+cHi6pv1C458k+UaSTw3+vzBjq5XkakZHQs8whbkN6n19UvNLckGSA20O\nT1fV85Oa2yK1JjKvZkVesj1bgmRS96B/s6reCtwKvC/J215TdHTcdqraXfu1jPF7PcToVYTrge8B\nfzHOwZO8Cfgs8P6q+vHws0nMrdX7x1bvJ0xoflX1s6q6ntE7X7+V5OZ5n49tbgvU2sqE5rWSL9me\nLUEy/2W/Dbw2Ec9IVX2v/fnfwD8xOlWZTbIeYBkvHC76YuEpnM74p/1C41BVvVwNo0PZuVOx7lpJ\nLmYUIo9W1WOTntug3t/N1Zvk/Nr4PwQ+z+jC4kR/b4NavzbBec29ZPsdYDfw28OXbCc1t7kJrvgP\nowtRLzG64LOGMVxsBd7I6DwP4BeAfwFu4QxeOFyiztW8/mLrpF5onF/r8kH7A8A/jKNW++wR4JPz\n+icyt1PUG/v8gLfQ7loweu/rq8DvTGJup6i1fhK/t3m1b2IKL9n+vN5Kh8hg4rcyulp/BNgxhvGu\naf+RDjC6pbij9a8FvszCt8I+0uofBn5vGTV2A/8F/B+jazzvOZPxOXmr7QjwV8us9V5GX77ngG8A\njzE6Fx5HrRsZnWMf4OQtym0TnNtC9W6dxPwYve/1H63Wc8CHzvTvRUetifze5tW+iZN3bSbyexv+\n+Ii8pG5nyzUSSauYQSKpm0EiqZtBIqmbQSKpm0EiqZtBIqmbQSKp2/8DGvC7Q10t8Z0AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1138694d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(image.getArray(),norm=LogNorm());"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       ..., \n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mask.getArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       ..., \n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.]], dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance.getArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Printout the WCS and FITS Header metaData keys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('NAXIS',\n",
       " 'EQUINOX',\n",
       " 'RADESYS',\n",
       " 'CRPIX1',\n",
       " 'CRPIX2',\n",
       " 'CD1_1',\n",
       " 'CD1_2',\n",
       " 'CD2_1',\n",
       " 'CD2_2',\n",
       " 'CRVAL1',\n",
       " 'CRVAL2',\n",
       " 'CUNIT1',\n",
       " 'CUNIT2',\n",
       " 'CTYPE1',\n",
       " 'CTYPE2')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metaWCS = wcs.getFitsMetadata()\n",
    "metaWCS.getOrderedNames()        #print them out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are  322 items.\n",
      "The first 20 names are ('VERSION', 'CREATOR', 'QEVAR', 'OVRDEP', 'SF', 'NB', 'NBULK', 'CHPANG', 'MAXY', 'MINY', 'MAXX', 'MINX', 'FPFILE', 'NF', 'PAIRID', 'PIXSIZE', 'DOMESEE', 'LASCPR', 'AERIND', 'H2OPRESS')\n"
     ]
    }
   ],
   "source": [
    "print \"There are \", metaData.nameCount(), \"items.\"\n",
    "print \"The first 20 names are\", metaData.names()[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's retrieve one of the header members and see who made this FITS file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PHOSIM\n"
     ]
    }
   ],
   "source": [
    "print metaData.get('CREATOR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics without a mask\n",
    "Now let's look at statistics of the image.  In a bit we will see what using a mask can do.  First set up the statistics flag:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The statistics flags are set to 0b110000001111.\n",
      "Errors will be calculated.\n"
     ]
    }
   ],
   "source": [
    "statFlags = math.NPOINT | math.MEAN | math.STDEV | math.MAX | math.MIN | math.ERRORS\n",
    "print \"The statistics flags are set to %s.\"%bin(statFlags)\n",
    "print \"Errors will be calculated.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's set individual pixel (0,0) to 65,000.  This will saturate a single pixel and screw up our overall statitics.  Later we will see how to mask it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "image.set(0, 0, 65000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*VERY IMPORTANT*: If you want to use a mask, you need to make a control object and tell the statistics routines which planes to pay attention to.  Also: unlike what is says in some documentation 0 will *ignore* all of the mask layers not use all of them. Here I will tell AFW to use the saturation (SAT) plane. Note I have't actually set any saturated bits yet.  I will do that later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "control = math.StatisticsControl()\n",
    "SAT = afwImg.MaskU_getPlaneBitMask(\"SAT\")\n",
    "control.setAndMask(SAT);        #pixels with this mask bit set will be ignored."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's do the statistics on the sensor.  The saturated pixel will screw things up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The image has dimensions 4000 x 4072 pixels\n",
      "Number of analyzed bins in image is 16288000\n",
      "Max    =     65000\n",
      "Min    =         0\n",
      "Mean   = 0.00555734 +- 0.0\n",
      "StdDev =     16.13\n"
     ]
    }
   ],
   "source": [
    "imageStatistics = math.makeStatistics(maskedImage, statFlags, control)\n",
    "numBins         = imageStatistics.getResult(math.NPOINT)[0]\n",
    "mean            = imageStatistics.getResult(math.MEAN)[0]\n",
    "\n",
    "print \"The image has dimensions %i x %i pixels\" %(maskedImage.getWidth(), maskedImage.getHeight())\n",
    "print \"Number of analyzed bins in image is %i\"  %numBins\n",
    "print \"Max    = %9d\"                            %imageStatistics.getResult(math.MAX)[0]\n",
    "print \"Min    = %9d\"                            %imageStatistics.getResult(math.MIN)[0]\n",
    "print \"Mean   = %9.8f +- %3.1f\"                 %imageStatistics.getResult(math.MEAN)\n",
    "print \"StdDev = %9.2f\"                          %imageStatistics.getResult(math.STDEV)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics with a mask\n",
    "Now also set the mask bit to SAT (Saturated) for that pixel. This should cause the builtin statistics routines to igrore this pixel as earlier we told them to pay attention to the SAT mask layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "saturationBit = mask.getPlaneBitMask('SAT')\n",
    "mask.set (0, 0, saturationBit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the let's try the same thing again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The image has dimensions 4000 x 4072 pixels\n",
      "Number of analyzed bins in image is 16287999\n",
      "Max    =      1103\n",
      "Min    =         0\n",
      "Mean   = 0.00156667 +- 0.0\n",
      "StdDev =      0.87\n"
     ]
    }
   ],
   "source": [
    "imageStatistics = math.makeStatistics(maskedImage, statFlags, control)\n",
    "numBins         = imageStatistics.getResult(math.NPOINT)[0]\n",
    "mean            = imageStatistics.getResult(math.MEAN)[0]\n",
    "\n",
    "print \"The image has dimensions %i x %i pixels\" %(maskedImage.getWidth(), maskedImage.getHeight())\n",
    "print \"Number of analyzed bins in image is %i\"  %numBins\n",
    "print \"Max    = %9d\"                            %imageStatistics.getResult(math.MAX)[0]\n",
    "print \"Min    = %9d\"                            %imageStatistics.getResult(math.MIN)[0]\n",
    "print \"Mean   = %9.8f +- %3.1f\"                 %imageStatistics.getResult(math.MEAN)\n",
    "print \"StdDev = %9.2f\"                          %imageStatistics.getResult(math.STDEV)[0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see one less pixel was considered and now the saturated pixel was ignored."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source Detection\n",
    "Now let's try to find the peak.  No backgound subtraction will be done because there is no background in this simulated image. This first step finds a \"footprint\" of connected pixels including its center. \n",
    "\n",
    "I am going to start by running these algorithms for detection and measurment \"manually\" so we can see the proper way to call them.  Later I will call routines that will call them for me."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "thresholdValue = 50\n",
    "npixMin        = 5 \n",
    "grow           = 1\n",
    "isotropic      = False\n",
    "\n",
    "threshold = afwDetect.Threshold(thresholdValue, afwDetect.Threshold.VALUE)\n",
    "footPrintSet  = afwDetect.FootprintSet(maskedImage, threshold, \"DETECTED\", npixMin)\n",
    "footPrintSet  = afwDetect.FootprintSet(footPrintSet, grow, isotropic)\n",
    "\n",
    "footPrints = footPrintSet.getFootprints()\n",
    "\n",
    "footPrintSet.setMask(maskedImage.getMask(), \"DETECTED\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's print out the peak we found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I found  1 footPrint(s)\n",
      "\n",
      "Footprint: 0\n",
      "A peak of value 1103.0\n",
      "was found at X = 2000.0 Y = 2032.0\n"
     ]
    }
   ],
   "source": [
    "print \"I found \", footPrints.size(), \"footPrint(s)\"\n",
    "print\n",
    "\n",
    "for i in range(0,footPrints.size()):\n",
    "    print \"Footprint:\",i\n",
    "    \n",
    "    peak = footPrints[i].getPeaks()[0]\n",
    "    print \"A peak of value\", peak.getPeakValue()\n",
    "    print \"was found at X =\", peak.getFx(),\"Y =\",peak.getFy()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We found our star!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source Measurement\n",
    "Now let's take that peak and try to measure it with one of the alogorithms.  We will use some of the infrastructure and also reduce the list of algorithms to apply to only one.  \n",
    "\n",
    "To do this I will also reset some \"aliases\" so that they don't point to algorithms that we aren't going to use (the default code uses several).  The aliases are pointers into records in the output table. The names of these pointers are independent of the algorithms.  So I can point them at the alogorithm results I want and then just ask for an output variable through an access function and it will get the correct one. This way, you can switch which algorithm generates the output, but not change how you access the answer.\n",
    "\n",
    "Note: This example uses tables and catalogs.  If you don't understand those you should also look at the short notebook which introduces them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "(1999.5, 2032.5)\n"
     ]
    }
   ],
   "source": [
    "# Make a Schema which is a list of columns for a table.  The object below store things in a table. \n",
    "# The Schema is a describption of what is in each column. It can be added to by algorithms later.\n",
    "schema = afwTable.SourceTable.makeMinimalSchema()\n",
    "\n",
    "# Make an object which we are going to use to configure the algorithms\n",
    "# Decide which algorithms to include (meas_algorithms/../tests/measure.py) is a good reference\n",
    "# Note: the pre-centroid fit is set seperately (it is set by default to 'centroid.sdss')\n",
    "measureSourcesConfig                  = measAlg.SourceMeasurementConfig()\n",
    "\n",
    "measureSourcesConfig.algorithms.names = ['shape.sdss']\n",
    "\n",
    "#Aliases only\n",
    "measureSourcesConfig.slots.psfFlux    = None\n",
    "measureSourcesConfig.slots.apFlux     = None\n",
    "measureSourcesConfig.slots.modelFlux  = None\n",
    "measureSourcesConfig.slots.instFlux   = None\n",
    "measureSourcesConfig.validate()\n",
    "\n",
    "# Now make the object which includes list the algorithms we want to use\n",
    "# using the measureSourcesConfig object I made before.\n",
    "# It also adds output places for those algorithms to the schema.\n",
    "# This object can apply it's algorithms to images given a list of sources.\n",
    "measureSources = measureSourcesConfig.makeMeasureSources(schema)\n",
    "\n",
    "# First make a catalog. Then, setup the aliases in the table so that it agrees with \n",
    "# our rules above (where we set the aliases to None).\n",
    "catalog  = afwTable.SourceCatalog(schema)\n",
    "measureSourcesConfig.slots.setupTable(catalog.getTable())\n",
    "\n",
    "# Take the set of footPrints (areas around the detected objects) we found when \n",
    "# we did the detection and then put them in the catalog.\n",
    "footPrintSet.makeSources(catalog)\n",
    "\n",
    "# Loop over all of the sources in the catalog.  For each source apply our measurement algorithms to it.\n",
    "# It uses the data in the exposure for the calculation.  After the algorithim is run, print the centroid\n",
    "# of the source as found by the algorithm.  We access it throgh an alias.\n",
    "for i, source in enumerate(catalog):\n",
    "    print i\n",
    "    measureSources.apply(source, exposure)\n",
    "    print source.getCentroid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source Detection and Measurement (the effecient coding way)\n",
    "Now let's do both steps again (detection and measurement) in a simpler way.  We will use the tasks that were setup to help us do these jobs.  This is the recomended interface. Internally it is doing a lot of what we did in the previous examples but being more comprehensive. We can configure them how we want.  Doing it this way is more compact, flexible and it is easier to understand since the steps pretty much make sense. BTW, you can still restrict the algorithms as I showed above etc in the same way. Here I use all of the default algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1999.5, 2032.5)\n",
      "     centroid.sdss (1999.5, 2032.5)\n",
      "     shape.sdss (ixx=2.88751926897, iyy=3.05814257946, ixy=-0.133327694904)\n",
      "     shape.sdss.centroid (1999.5, 2032.5)\n",
      "     flux.gaussian 21410.8483803\n"
     ]
    }
   ],
   "source": [
    "# Configure the detection and measurement algorithms\n",
    "schema                = afwTable.SourceTable.makeMinimalSchema()\n",
    "detectSourcesConfig   = measAlg.SourceDetectionConfig(thresholdType='value')\n",
    "measureSourcesConfig  = measAlg.SourceMeasurementConfig()\n",
    "\n",
    "# Setup the detection and measurement tasks\n",
    "detect  = measAlg.SourceDetectionTask  (config=detectSourcesConfig,  schema=schema)\n",
    "measure = measAlg.SourceMeasurementTask(config=measureSourcesConfig, schema=schema)\n",
    "\n",
    "# Set flux aliases to None; a hack for an incompatability between makeMinimalSchema() and the\n",
    "# default SourceMeasurementConfig() options.\n",
    "measureSourcesConfig.slots.psfFlux    = None\n",
    "measureSourcesConfig.slots.apFlux     = None\n",
    "measureSourcesConfig.slots.modelFlux  = None\n",
    "measureSourcesConfig.slots.instFlux   = None\n",
    "measureSourcesConfig.validate()\n",
    "\n",
    "# Detect the sources,then put them into a catalog (the table is where the catalog atually stores stuff)\n",
    "table   = afwTable.SourceTable.make(schema)\n",
    "catalog = detect.makeSourceCatalog(table, exposure, sigma=5)\n",
    "\n",
    "# Get the sources out of the catalog\n",
    "sources = catalog.sources\n",
    "\n",
    "# Apply the measurement routines to the exposure using the sources as input\n",
    "measure.run(exposure, sources)\n",
    "\n",
    "# Now let's look at the output from some of the measurment algorithms.\n",
    "fields = ['centroid.sdss', 'shape.sdss', 'shape.sdss.centroid','flux.gaussian']\n",
    "keys   = [schema.find(f).key for f in fields]\n",
    "\n",
    "for source in sources:\n",
    "    print source.getCentroid() #This uses one of the aliases\n",
    "\n",
    "    # Now loop through the keys we want\n",
    "    for f,k in zip(fields, keys):\n",
    "        print '    ', f, source.get(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the DS9 program to display our results.  This will overlay the footprint and the source measured centroid to the sensor data. This is commented out by default in case you don't have ds9 installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#import lsst.afw.display.ds9 as ds9\n",
    "##ds9.setMaskPlaneVisibility(\"DETECTED\",False)\n",
    "#ds9.setMaskTransparency(75)\n",
    "\n",
    "#ds9.mtv(maskedImage)\n",
    "#ds9.dot(\"+\", source.getX(), source.getY())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can take the code in this tutorial and turn it into a single program.  Then, we can learn how to make it into a commandLine Task and use some framework code that will take care of some more things I did in the last step for us automatically.  I do that in the next notebook."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}