photometry.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "from astropy.io import fits\n",
    "from astropy.stats import sigma_clip\n",
    "from astropy.coordinates import SkyCoord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# reset defalult plotting values\n",
    "plt.rcParams['figure.figsize'] = (10, 7)\n",
    "plt.rc('font', family='sans-serif')\n",
    "plt.rc('axes', labelsize=14)\n",
    "plt.rc('axes', labelweight='bold')\n",
    "plt.rc('axes', titlesize=16)\n",
    "plt.rc('axes', titleweight='bold')\n",
    "plt.rc('axes', linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Photometry\n",
    "## quantifying light from astronomical sources\n",
    "\n",
    "<img src=\"media/photometry.png\" width=400>\n",
    "\n",
    "### Prof. Robert Quimby\n",
    "&copy; 2019 Robert Quimby"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## In this tutorial you will..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Learn what it means to do aperture photometry on a celestial target\n",
    "- measure the number of object counts above the background level recorded by a CCD\n",
    "- determine the uncertainty in this measurement\n",
    "- use the `astropy` package, `photutils` to measure counts in a circular aperture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Terminology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](http://i-fiberoptics.com/images/if-pm.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the context of photometry, \"aperture\" refers **not** to the size of your telescope, but to the area of the telescope's focal plane used when measuring light."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simplest form: box aperture photometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load (processed) image data\n",
    "image = fits.getdata('media/example.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# show a close-up view of a star\n",
    "x, y = 244, 640\n",
    "margin = 10\n",
    "\n",
    "# show the stamp\n",
    "stamp = image[y - margin : y + margin + 1, x - margin : x + margin + 1]\n",
    "xorigin, yorigin = x - margin, y - margin # for later use\n",
    "extent = [x - margin,  x + margin + 1, y - margin, y + margin + 1]\n",
    "plt.imshow(stamp, origin='lower', cmap='gray', vmin=1350, vmax=1700, extent=extent);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# define a background aperture\n",
    "bgx, bgy = x, y + 20\n",
    "bgstamp = image[bgy - margin : bgy + margin + 1, bgx - margin : bgx + margin + 1]\n",
    "bgextent = [bgx - margin,  bgx + margin + 1, bgy - margin, bgy + margin + 1]\n",
    "plt.imshow(bgstamp, origin='lower', cmap='gray', vmin=1350, vmax=1700, extent=bgextent);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# total counts in the target aperture\n",
    "counts = ????\n",
    "\n",
    "# total counts in the background aperture\n",
    "bgcounts = ????\n",
    "\n",
    "print(\"    target aperture has {:.0f} counts\".format(counts))\n",
    "print(\"background aperture has {:.0f} counts\".format(bgcounts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# how many counts do we get from just the target?\n",
    "print(\"That is {:.0f} counts from the object\".format(????))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What is the Uncertainty in the Object Counts?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recall Propagation of Uncertainties\n",
    "\n",
    "* Given a function $f$ of several variables, $x_0, x_1 ... x_{n-1}$, if the variables have random, independent uncertainties, $\\sigma_{x_0}, \\sigma_{x_1} ... \\sigma_{x_{n-1}}$, then the uncertainty in the function, $\\sigma_f$, is:\n",
    "\n",
    "$$ \\sigma_f = \\sqrt{ \\left( {\\partial f \\over \\partial x_0 }\\sigma_{x_0}\\right)^2\n",
    "                   + \\left( {\\partial f \\over \\partial x_1 }\\sigma_{x_1}\\right)^2\n",
    "                   + ...\n",
    "                   + \\left( {\\partial f \\over \\partial x_n }\\sigma_{x_{n-1}}\\right)^2\n",
    "                    } $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The object counts, $f$, are computed from the equation:\n",
    "$$f = \\sum_{i=0}^{N-1} (S_i + B_i - \\bar{B})$$\n",
    "where:\n",
    " * $N$ is the number of pixels (i.e. the area) in the aperture\n",
    " * $S_i$ is the source (star) counts in the $i^{\\rm th}$ pixel,\n",
    " * $B_i$ is the background counts in the $i^{\\rm th}$ pixel, \n",
    " * $\\bar{B}$ is the average value of the $M$ pixels in the background sample\n",
    " \n",
    "$$\\bar{B} = {1 \\over M}\\sum_{j=0}^{M-1} B_j$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Estimate the uncertainty in $\\bar{B}$ from the background sample.\n",
    "\n",
    "Assume that the total source counts, $S = \\sum S_i$, are drawn from a Poisson distribution. Then using the fact that Poisson distributions have $\\sigma^2 = \\mu$, we find:\n",
    "\n",
    "$$\\sigma^2_S = S$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So the uncertainty in the measured object counts, $\\sigma_f$ is given by:\n",
    "$$\\sigma_f = \\sqrt{ S + N\\sigma_B^2 + {N^2\\sigma_B^2 \\over M} } $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Calculate the uncertainty in the counts measurement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# photometric uncertainty\n",
    "obj_var = ????\n",
    "bg_var = stamp.size * bgstamp.var()\n",
    "avbg_var = bgstamp.std()**2 / bgstamp.size * stamp.size**2\n",
    "ecounts = np.sqrt( obj_var + bg_var +  avbg_var)\n",
    "print(\"Total uncertainty is {:.0f} counts\".format(ecounts))\n",
    "print(\"({:.0f} object, {:.0f} background, {:.0f} background average)\".format(obj_var**0.5, bg_var**0.5, avbg_var**0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Background pixels just add noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_example_star_image(showaperture=False):\n",
    "    star = fits.getdata('media/star_image.fits')\n",
    "    plt.figure(figsize=(10, 10))\n",
    "    plt.imshow(star, origin='lower')\n",
    "    ny, nx = star.shape\n",
    "    for x in range(nx):\n",
    "        for y in range(ny):\n",
    "            plt.annotate(\"{:.0f}\".format(star[y, x]), xy=(x, y), va='bottom', ha='center')\n",
    "            \n",
    "    if showaperture:\n",
    "        cx, cy = 5.2, 5.6\n",
    "        circle= plt.Circle((cx, cy), radius=3, fill=False, color='r', lw=2)\n",
    "        ax=plt.gca()\n",
    "        ax.add_patch(circle);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plot_example_star_image()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Circular Aperture Photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- sum up the light from a star in a circular aperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plot_example_star_image(showaperture=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Decide where to center the aperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# get the x,y coordinates for each pixel in the stamp\n",
    "yinds, xinds = np.indices(stamp.shape)\n",
    "\n",
    "# calculate the \"weight\" for each pixel\n",
    "weight = ????\n",
    "\n",
    "## calculate the \"center of mass\" in each coordinate\n",
    "cmx = ????\n",
    "cmy = ????\n",
    "print(\"object is centered at {:.2f}, {:.2f}\".format(cmx, cmy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# mark the intensity weighted center on the image stamp\n",
    "plt.imshow(stamp, origin='lower', cmap='gray', vmin=1350, vmax=1700)\n",
    "plt.plot(cmx, cmy, 'r+', ms=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Decide on the radius for the aperture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Light is spread out due to refracton and atomspheric seeing\n",
    "- cannot measure 100% of the light in a finite radius aperture\n",
    "- but you can measure as close to this as you care to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# plot the object's profile along one axis\n",
    "profile = ????\n",
    "profile -= ????\n",
    "profile /= ????\n",
    "plt.plot(profile)\n",
    "plt.axvline(cmx, color='r', ls='dashed')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(profile)\n",
    "plt.axvline(cmx, color='r', ls='dashed')\n",
    "fwhm = ????\n",
    "plt.plot([cmx - fwhm/2, cmx + fwhm/2], [0.5, 0.5], color='red', ls='dotted')\n",
    "\n",
    "plt.xlabel('Stamp Column Number')\n",
    "plt.ylabel('Column Sum');\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# show the circular aperture on the image stamp\n",
    "plt.imshow(stamp, origin='lower', cmap='gray', vmin=1350, vmax=5700)\n",
    "plt.plot(cmx, cmy, 'r+', ms=14);\n",
    "\n",
    "circle = plt.Circle((cmx, cmy), radius=fwhm, fill=False, color='r', lw=2)\n",
    "ax = plt.gca()\n",
    "ax.add_patch(circle);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sum up the Counts in a Circular Aperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plot_example_star_image(showaperture=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Rules for measuring counts in an aperture\n",
    "* if a pixel is wholly contained in the aperture, add all of its counts\n",
    "* if $XX$% of a pixel is inside the aperture, add $XX$% of the counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Use the `photutils` package for the dirty work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# sum the counts in the aperture\n",
    "import photutils\n",
    "aperture = photutils.CircularAperture( ???? )\n",
    "phot_table = photutils.aperture_photometry(????)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# look at the photometry table\n",
    "phot_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Define a background sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# show the circular aperture on the image stamp\n",
    "cmx2, cmy2 = ????\n",
    "\n",
    "plt.imshow(image, origin='lower', cmap='gray', vmin=1350, vmax=1700)\n",
    "plt.xlim(cmx2 - 7 * fwhm, cmx2 + 7 * fwhm)\n",
    "plt.ylim(cmy2 - 7 * fwhm, cmy2 + 7 * fwhm)\n",
    "\n",
    "plt.plot(cmx2, cmy2, 'r+', ms=14);\n",
    "circle = plt.Circle((cmx2, cmy2), radius=fwhm, fill=False, color='r', lw=2)\n",
    "ax = plt.gca()\n",
    "ax.add_patch(circle);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Get the background sample pixels in an annulus around the target position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# calculate the distance from each image pixel to the centroid\n",
    "imy, imx = ????\n",
    "dists = ????\n",
    "wbg = ????\n",
    "bg_sample = ????"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# measure the background subtracted counts in the aperture\n",
    "phot_table = photutils.aperture_photometry(????)\n",
    "print(\"{:.0f} counts from the object\".format(????))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Determine the Uncetainty in the Circular Aperture Counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "obj_var = ????\n",
    "bg_var = ???? * bg_sample.var()\n",
    "avbg_var = bg_sample.var() / bg_sample.size * ????**2\n",
    "ecounts = np.sqrt(obj_var + bg_var +  avbg_var)\n",
    "print(\"Total uncertainty is {:.0f} counts\".format(ecounts))\n",
    "print(\"({:.0f} object, {:.0f} background, {:.0f} background average)\".format(obj_var**0.5, bg_var**0.5, avbg_var**0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More on photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[`photutils.aperture` photometry documentation](https://photutils.readthedocs.io/en/stable/aperture.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.4"
  },
  "livereveal": {
   "scroll": "True",
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}