03_aperture_photometry.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"title_ID\"></a>\n",
    "\n",
    "<a href=\"http://photutils.readthedocs.io/en/stable/index.html\"><img src=\"https://photutils.readthedocs.io/en/stable/_static/photutils_banner.svg\" width=300></a>\n",
    "\n",
    "# Aperture Photometry with `photutils`\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### What is aperture photometry?\n",
    "The most common method to measure the flux of a celestial source is aperture photometry. This kind of photometry measures the amount of flux within a region of an astronomical image of defined shape and size (an aperture) surrounding a source. The ideal aperture would capture all of the flux emitted by a desired source, and none of the flux emitted by the surrounding sky or nearby sources. Especially when performing photometry on image data that includes a number of sources with varying size and shape, it is important to perform aperture corrections to account for imperfect apertures and better constrain photometric errors.\n",
    "\n",
    "The `photutils` package provides tools for performing photometry with apertures of various shapes.\n",
    "\n",
    "##### What does this tutorial include?\n",
    "This tutorial covers how to perform aperture photometry with `photutils`, including the following methods:\n",
    "* Creating Apertures\n",
    "    * Circular Apertures\n",
    "    * Elliptical Apertures\n",
    "    * Sky Apertures with WCS\n",
    "* Performing Aperture Photometry\n",
    "* Calculating Aperture Corrections with Local Background Subtraction\n",
    "\n",
    "##### Which data are used in this tutorial?\n",
    "We will be manipulating Hubble eXtreme Deep Field (XDF) data, which was collected using the Advanced Camera for Surveys (ACS) on Hubble between 2002 and 2012. The image we use here is the result of 1.8 million seconds (500 hours!) of exposure time, and includes some of the faintest and most distant galaxies that have ever been observed. \n",
    "\n",
    "*The methods demonstrated here are available in narrative form within the `photutils.aperture` [documentation]( http://photutils.readthedocs.io/en/stable/aperture.html).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<b>Important:</b> Before proceeding, please be sure to update your versions of <code>astropy</code>, <code>matplotlib</code>, and <code>photutils</code>, or this notebook may not work properly. Or, if you don't want to handle packages individually, you can always use (and keep updated!) the <a href=\"https://astroconda.readthedocs.io\">AstroConda</a> distribution.\n",
    " \n",
    "</div>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import necessary packages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import packages that we will use to perform arithmetic functions and visualize data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "import astropy.units as u\n",
    "from astropy.nddata import CCDData\n",
    "from astropy.stats import sigma_clipped_stats, SigmaClip\n",
    "from astropy.visualization import ImageNormalize, LogStretch\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import LogLocator\n",
    "import numpy as np\n",
    "from photutils.background import Background2D, MeanBackground\n",
    "\n",
    "# Show plots in the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also define some `matplotlib` parameters, such as title font size and the dpi, to make sure our plots look nice. To make it quick, we'll do this by loading a [style file shared with the other photutils tutorials](../photutils_notebook_style.mplstyle) into `pyplot`. We will use this style file for all the notebook tutorials. (See [here](https://matplotlib.org/users/customizing.html) to learn more about customizing `matplotlib`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('../photutils_notebook_style.mplstyle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieve data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the introduction, we will be using Hubble eXtreme Deep Field (XDF) data. Since this file is too large to store on GitHub, we will just use `astropy` to directly download the file from the STScI archive: https://archive.stsci.edu/prepds/xdf/ \n",
    "\n",
    "(Generally, the best package for web queries of astronomical data is [Astroquery](https://astroquery.readthedocs.io/en/latest/); however, the dataset we are using is a High Level Science Product (HLSP) and thus is not located within a catalog that could be queried with Astroquery.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'https://archive.stsci.edu/pub/hlsp/xdf/hlsp_xdf_hst_acswfc-60mas_hudf_f435w_v1_sci.fits'\n",
    "with fits.open(url) as hdulist:\n",
    "    hdulist.info()\n",
    "    data = hdulist[0].data\n",
    "    header = hdulist[0].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As explained in a [previous notebook](../01_background_estimation/01_background_estimation.ipynb) on background estimation, it is important to **mask** these data, as a large portion of the values are equal to zero. We will mask out the non-data portions of the image array, so all of those pixels that have a value of zero don't interfere with our statistics and analyses of the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the mask\n",
    "mask = data == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Throughout this notebook, we are going to store our images in Python using a `CCDData` object (see [Astropy documentation](http://docs.astropy.org/en/stable/nddata/index.html#ccddata-class-for-images)), which contains a `numpy` array in addition to metadata such as uncertainty, masks, or units. In this case, our data is in electrons (counts) per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = u.electron / u.s\n",
    "xdf_image = CCDData(data, unit=unit, meta=header, mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Set up the figure with subplots\n",
    "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
    "\n",
    "# Set up the normalization and colormap\n",
    "norm_image = ImageNormalize(vmin=1e-4, vmax=5e-2, stretch=LogStretch(), clip=False)\n",
    "cmap = plt.get_cmap('viridis')\n",
    "cmap.set_over(cmap.colors[-1])\n",
    "cmap.set_under(cmap.colors[0])\n",
    "cmap.set_bad('white') # Show masked data as white\n",
    "xdf_image_clipped = np.clip(xdf_image, 1e-4, None) # clip to plot with logarithmic stretch\n",
    "\n",
    "# Plot the data\n",
    "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n",
    "                      norm=norm_image, cmap=cmap)\n",
    "\n",
    "# Define the colorbar and fix the labels\n",
    "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
    "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
    "cbar.ax.set_yticklabels(labels)\n",
    "\n",
    "# Define labels\n",
    "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n",
    "               rotation=270, labelpad=30)\n",
    "ax1.set_xlabel('X (pixels)')\n",
    "ax1.set_ylabel('Y (pixels)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Tip: Double-click on any inline plot to zoom in.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Creating Apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `photutils`, users can create apertures with the following shapes:\n",
    "\n",
    "<a><img src=\"apertures.png\" width=700 alt=\"Examples of circular, elliptical, and rectangular apertures and annuli.\"></a>\n",
    "\n",
    "Each of these can be defined either in pixel coordinates or in celestial coordinates (using a WCS transformation).\n",
    "\n",
    "It is also possible for users to create custom aperture shapes.\n",
    "\n",
    "Any aperture object is created by defining its position and size (and, if applicable, orientation). Let's use the `find_peaks` method that we learned about in a [previous notebook](../02_source_detection/02_source_detection.ipynb) to get the positions of sources in our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils import find_peaks\n",
    "from photutils.centroids import centroid_2dg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate statistics\n",
    "mean, median, std = sigma_clipped_stats(xdf_image.data, sigma=3.0, maxiters=5, mask=xdf_image.mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources_findpeaks = find_peaks(xdf_image.data, mask=xdf_image.mask, \n",
    "                               threshold=20.*std, box_size=30, \n",
    "                               centroid_func=centroid_2dg)     \n",
    "# Display the table\n",
    "sources_findpeaks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's plot the centroids of each of these sources:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the figure with subplots\n",
    "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
    "\n",
    "# Plot the data\n",
    "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), norm=norm_image)\n",
    "ax1.scatter(sources_findpeaks['x_centroid'], sources_findpeaks['y_centroid'], s=10, marker='.', \n",
    "            lw=1, alpha=0.7, color='r')#facecolor='None', edgecolor='r')\n",
    "\n",
    "# Define the colorbar\n",
    "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
    "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
    "cbar.ax.set_yticklabels(labels)\n",
    "\n",
    "# Define labels\n",
    "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n",
    "               rotation=270, labelpad=30)\n",
    "ax1.set_xlabel('X (pixels)')\n",
    "ax1.set_ylabel('Y (pixels)')\n",
    "ax1.set_title('find\\_peaks Sources')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So thanks to `find_peaks`, we now we know the positions of all our sources. Next, we need to define apertures for each source. First, as the simplest example, let's try using circular apertures of a fixed size.\n",
    "\n",
    "### Circular Apertures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils import CircularAperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the aperture\n",
    "position = (sources_findpeaks['x_centroid'], sources_findpeaks['y_centroid'])\n",
    "radius = 10.\n",
    "circular_aperture = CircularAperture(position, r=radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the figure with subplots\n",
    "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
    "\n",
    "# Plot the data\n",
    "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), norm=norm_image)\n",
    "\n",
    "# Plot the apertures\n",
    "circular_aperture.plot(color='red', alpha=0.7)\n",
    "\n",
    "# Define the colorbar\n",
    "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
    "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
    "cbar.ax.set_yticklabels(labels)\n",
    "\n",
    "# Define labels\n",
    "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n",
    "               rotation=270, labelpad=30)\n",
    "ax1.set_xlabel('X (pixels)')\n",
    "ax1.set_ylabel('Y (pixels)')\n",
    "ax1.set_title('Circular Apertures')\n",
    "\n",
    "# Crop to show an inset of the data\n",
    "ax1.set_xlim(2000, 3000)\n",
    "ax1.set_ylim(2000, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, these circular apertures don't fit our data very well. After all, this is the Hubble eXtreme Deep Field, so there aren't any nice, round, nearby Milky Way stars in this image! \n",
    "\n",
    "Let's use ellipses instead, to better match the morphology of the galactic blobs in our image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elliptical Apertures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils import (detect_sources, source_properties, \\\n",
    "                       EllipticalAnnulus, EllipticalAperture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a [previous notebook](../02_source_detection/02_source_detection.ipynb), we showed how you can use the `photutils.detect_sources` [feature](https://photutils.readthedocs.io/en/stable/api/photutils.detect_sources.html) to generate segmentation maps, which identify and label contiguous (connected) objects within an image. Then, with `source_properties` [feature](https://photutils.readthedocs.io/en/stable/api/photutils.segmentation.source_properties.html?highlight=source_properties), you can access descriptive properties for each unique object - not just their centroid positions, but also their pixel areas, eccentricities, orientations with respect to the coordinate frame of the image, and more.\n",
    "\n",
    "Here we'll use the centroid, semimajor axis, semiminor axis, and orientation values from `source_properties` to generate elliptical apertures for each of the sources in our image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Define threshold and minimum object size\n",
    "threshold = 5. * std\n",
    "npixels = 15\n",
    "\n",
    "# Create a segmentation image\n",
    "segm = detect_sources(xdf_image.data, threshold, npixels)\n",
    "\n",
    "# Create a catalog using source properties\n",
    "catalog = source_properties(xdf_image.data, segm)\n",
    "table = catalog.to_table()\n",
    "\n",
    "# Display the table\n",
    "table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 3.  # approximate isophotal extent of semimajor axis\n",
    "\n",
    "# Create the apertures\n",
    "elliptical_apertures = []\n",
    "for obj in catalog:\n",
    "    position = (obj.xcentroid.value, obj.ycentroid.value)\n",
    "    a = obj.semimajor_axis_sigma.value * r\n",
    "    b = obj.semiminor_axis_sigma.value * r\n",
    "    theta = obj.orientation.value\n",
    "    elliptical_apertures.append(EllipticalAperture(position, a, b, theta=theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the figure with subplots\n",
    "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
    "\n",
    "# Plot the data\n",
    "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), norm=norm_image)\n",
    "\n",
    "# Plot the apertures\n",
    "for aperture in elliptical_apertures:\n",
    "    aperture.plot(color='red', alpha=0.7, ax=ax1)\n",
    "\n",
    "# Define the colorbar\n",
    "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
    "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
    "cbar.ax.set_yticklabels(labels)\n",
    "\n",
    "# Define labels\n",
    "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n",
    "               rotation=270, labelpad=30)\n",
    "ax1.set_xlabel('X (pixels)')\n",
    "ax1.set_ylabel('Y (pixels)')\n",
    "ax1.set_title('Elliptical Apertures')\n",
    "\n",
    "# Crop to show an inset of the data\n",
    "ax1.set_xlim(2000, 3000)\n",
    "ax1.set_ylim(2000, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, these custom-made elliptical apertures fit our XDF galaxies much better than the one-size-fits-all circular apertures from before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sky Coordinates & Apertures\n",
    "\n",
    "At the moment, the positions of our apertures are in pixels, relative to our data array. However, if you need aperture positions in terms of celestial coordinates, `photutils` also includes aperture objects that can be integrated with Astropy's `SkyCoords`.\n",
    "\n",
    "Fortunately this is extremely easy when we use the [World Coordinate System (WCS)](http://docs.astropy.org/en/stable/wcs/) to produce a WCS object from the header of the FITS file containing our image, and then the `to_sky()` method to transform our `EllipticalAperture` objects into `SkyEllipticalAperture` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.wcs import WCS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wcs = WCS(header)\n",
    "sky_elliptical_apertures = [ap.to_sky(wcs) for ap in elliptical_apertures]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to generate `SkyEllipticalAperture` objects from the get-go, we could have used that WCS object in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils import SkyEllipticalAperture\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "r = 3.  # approximate isophotal extent of semimajor axis\n",
    "\n",
    "# Create the apertures\n",
    "sky_elliptical_apertures = []\n",
    "for obj in catalog:\n",
    "    # Convert the centroids into RA/Dec using WCS\n",
    "    ra, dec = wcs.all_pix2world(obj.xcentroid.value, obj.ycentroid.value, 0)\n",
    "    # Convert the positions to an Astropy SkyCoord object, with units!\n",
    "    sky_position = SkyCoord(ra, dec, unit=u.deg)\n",
    "    \n",
    "    # Define the elliptical parameters, now with units\n",
    "    a = obj.semimajor_axis_sigma.value * r * u.pix\n",
    "    b = obj.semiminor_axis_sigma.value * r * u.pix\n",
    "    theta = obj.orientation.value  * u.rad\n",
    "    \n",
    "    # Convert the theta from radians from X axis to the radians from North \n",
    "    x_to_north_angle = (90. + header['ORIENTAT']) * u.deg\n",
    "    x_to_north_angle_rad = x_to_north_angle.to_value(u.rad) * u.rad\n",
    "    theta -= x_to_north_angle_rad\n",
    "    \n",
    "    # Define the apertures\n",
    "    ap = SkyEllipticalAperture(sky_position, a, b, theta=theta)\n",
    "    sky_elliptical_apertures.append(ap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, you can't plot SkyApertures. However, you can use them just as easily to perform aperture photometry!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## Performing Aperture Photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have aperture objects that fit our data reasonably well, we can finally perform photometry with the `aperture_photometry` function. This function takes the following arguments:\n",
    "\n",
    "* **`data`** - the background-subtracted data array on which to perform photometry.\n",
    "* **`apertures`** - an aperture object containing the aperture(s) to use for the photometry.\n",
    "* **`error`** (optional) - an array of values that represent the pixel-wise Gaussian 1-sigma errors of the input data.\n",
    "* **`mask`** (optional) - a mask for the `data` to exclude certain pixels from calculations.\n",
    "* **`method`** (optional) - how to place the aperture(s) onto the pixel grid (see below).\n",
    "* **`unit`** (optional) - unit of `data` and `error`.\n",
    "* **`wcs`** (optional) - the WCS transformation to use if `apertures` is a `SkyAperture` object. \n",
    "* **`subpixels`** (optional) - the factor by which pixels are resampled (see below).\n",
    "\n",
    "The following methods are the options for how to place apertures onto the data pixel grid:\n",
    "\n",
    "* **exact** (default) - calculate the exact fractional overlap of each aperture for each overlapping pixel. This method is the most accurate, but will also take the longest. \n",
    "* **center** - a pixel is either entirely in or entirely out of the aperture, depending on whether the pixel center is inside or outside of the aperture.\n",
    "* **subpixel** - a pixel is divided into `subpixels` x `subpixels` subpixels, each of which are considered to be entirely in or out of the aperture depending on whether its center is in or out of the aperture. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils import aperture_photometry\n",
    "from astropy.table import QTable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what this looks like using the first aperture in our image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The CCDData mask will be automatically applied\n",
    "phot_datum = aperture_photometry(xdf_image, elliptical_apertures[0]) \n",
    "phot_datum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `aperture_sum` value is what reports the number of electron counts within the aperture: 3.47 e<sup>–</sup>/s.\n",
    "\n",
    "And, just as a check, to make sure our sky apertures give basically the same answer..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The CCDData mask will be automatically applied\n",
    "sky_phot_datum = aperture_photometry(xdf_image, sky_elliptical_apertures[0], wcs=wcs)\n",
    "sky_phot_datum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Woohoo!\n",
    "\n",
    "Unfortunately for our purposes, the `aperture_photometry` function can be only used alone for one of the two cases:\n",
    "* Identical apertures at distinct positions (e.g. circular apertures with `r = 3` for many sources)\n",
    "* Distinct apertures at identical positions (e.g. two circular apertures with `r = 3` and `r = 5` for one source)\n",
    "\n",
    "Since our elliptical apertures are distinct apertures at distinct positions, we need to do a little more work to get a single table of photometric values.\n",
    "\n",
    "(This step will take a while, almost 5 minutes, so hang tight!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The CCDData mask will be automatically applied\n",
    "phot_table = aperture_photometry(xdf_image, elliptical_apertures[0])\n",
    "id = 1\n",
    "for aperture in elliptical_apertures[1:]:\n",
    "    id += 1\n",
    "    phot_row = aperture_photometry(xdf_image, aperture)[0]\n",
    "    phot_row[0] = id\n",
    "    phot_table.add_row(phot_row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at all these apertures we've made:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the table\n",
    "phot_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's only so much you can learn from looking at a table of numbers, so let's explore alternate ways to examine these data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 5))\n",
    "\n",
    "values = [phot.value for phot in phot_table['aperture_sum']]\n",
    "logbins=bins = 10.**(np.linspace(-1, 2, 100))\n",
    "plt.hist(values, bins=logbins)\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.title('Histogram of Source Photometry')\n",
    "plt.xlabel(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')))\n",
    "plt.ylabel('Number of Sources')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 5))\n",
    "\n",
    "plt.scatter(table['area'], values, alpha=0.5)\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.title('Count Rate v. Aperture Area')\n",
    "plt.xlabel('Aperture Area [pixels$^2$]')\n",
    "plt.xlabel(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "<h3>Exercise:</h3><br>\n",
    "\n",
    "Re-calculate the photometry for these elliptical apertures - or just a subset of them - using the <code>subpixel</code> aperture placement method instead of the default <code>exact</code> method. How does this affect the count sum calculated for those apertures?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Aperture Corrections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've done photometry with some lovely apertures, but unfortunately even using elliptical apertures with unique sizes and orientations does not account for an important source of extraneous flux: the sky background.\n",
    "\n",
    "## Local Background Subtraction\n",
    "\n",
    "In the [background estimation notebook](../01_background_estimation/01_background_estimation.ipynb), we explored how to perform global background subtraction of image data with `photutils`. However, you can also use `photutils` to perform local background estimations for aperture corrections.\n",
    "\n",
    "To estimate the local background for each aperture, measure the counts within annulus apertures around (but not including!) each source. In our example, we defined elliptical apertures with `r = 3` to measure the counts within each source. To calculate the background for each source, let's measure the counts elliptical annuli between `r = 3.5` and `r = 5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_in = 3.5  # approximate isophotal extent of inner semimajor axis\n",
    "r_out = 5.  # approximate isophotal extent of inner semimajor axis\n",
    "\n",
    "# Create the apertures\n",
    "elliptical_annuli = []\n",
    "for obj in catalog:\n",
    "    position = (obj.xcentroid.value, obj.ycentroid.value)\n",
    "    a_in = obj.semimajor_axis_sigma.value * r_in\n",
    "    a_out = obj.semimajor_axis_sigma.value * r_out\n",
    "    b_out = obj.semiminor_axis_sigma.value * r_out\n",
    "    theta = obj.orientation.value\n",
    "    elliptical_annuli.append(EllipticalAnnulus(position, a_in, a_out, b_out, theta=theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the figure with subplots\n",
    "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
    "\n",
    "# Plot the data\n",
    "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n",
    "                      norm=norm_image)\n",
    "\n",
    "# Plot the apertures\n",
    "for aperture in elliptical_annuli:\n",
    "    aperture.plot(color='red', alpha=0.4, ax=ax1, fill=True)\n",
    "for aperture in elliptical_apertures:\n",
    "    aperture.plot(color='white', alpha=0.7, ax=ax1)\n",
    "\n",
    "# Define the colorbar\n",
    "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
    "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
    "cbar.ax.set_yticklabels(labels)\n",
    "\n",
    "# Define labels\n",
    "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n",
    "               rotation=270, labelpad=30)\n",
    "ax1.set_xlabel('X (pixels)')\n",
    "ax1.set_ylabel('Y (pixels)')\n",
    "ax1.set_title('Elliptical Annuli')\n",
    "\n",
    "# Crop to show an inset of the data\n",
    "ax1.set_xlim(2000, 3000)\n",
    "ax1.set_ylim(2000, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating Aperture Corrections\n",
    "\n",
    "Now that our apertures have been defined, we can do photometry with them to estimate and account for the background. The aperture correction is calculated by:\n",
    "- Calculating the count rate within each annulus using `aperture_photometry`\n",
    "- Dividing each annulus' count rate by each annulus' area to get the mean background value for each annulus\n",
    "- Taking the mean of those annulus means to get a mean background value for the entire image\n",
    "- Multiplying the global background mean value times the area of each elliptical photometric aperture, to get the estimated background count rate within each aperture\n",
    "- Subtracting the estimated background count rate from the photometric count rate for each aperture\n",
    "\n",
    "(Just like when we did photometry with the elliptical apertures above, the below step will take almost 5 minutes.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The CCDData mask will be automatically applied\n",
    "bkg_phot_table = aperture_photometry(xdf_image, elliptical_annuli[0])\n",
    "id = 1\n",
    "for aperture in elliptical_annuli[1:]:\n",
    "    id += 1\n",
    "    phot_row = aperture_photometry(xdf_image, aperture)[0]\n",
    "    phot_row[0] = id\n",
    "    bkg_phot_table.add_row(phot_row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display table\n",
    "bkg_phot_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the mean background level (per pixel) in the annuli \n",
    "bkg_area = [annulus.area() for annulus in elliptical_annuli]\n",
    "bkg_mean_per_aperture = bkg_phot_table['aperture_sum'].value / bkg_area\n",
    "bkg_mean = np.average(bkg_mean_per_aperture) * (u.electron / u.s)\n",
    "print('Background mean:', bkg_mean)\n",
    "\n",
    "# Calculate the total background within each elliptical aperture\n",
    "bkg_sum = bkg_mean * table['area'].value\n",
    "\n",
    "# Subtract the background from the original photometry\n",
    "flux_bkgsub = phot_table['aperture_sum'] - bkg_sum\n",
    "\n",
    "# Add this as a column to the original photometry table\n",
    "phot_table['aperture_sum_bkgsub'] = flux_bkgsub"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might have noticed that these background count rates are *really* small. In this case, this is to be expected - since our example XDF data is a high-level science product (HLSP) that already has already been background-subtracted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's see the difference between our original count rates and our background-subtracted count rates (it should be small for us!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display table\n",
    "phot_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 5))\n",
    "\n",
    "values = [phot.value for phot in phot_table['aperture_sum']]\n",
    "values_bkgsub = [phot.value for phot in phot_table['aperture_sum_bkgsub']]\n",
    "logbins=bins = 10.**(np.linspace(-1, 2, 100))\n",
    "plt.hist(values, bins=logbins, alpha=0.7, label='Original photometry')\n",
    "plt.hist(values_bkgsub, bins=logbins, alpha=0.7, label='Background-subtracted')\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.title('Histogram of Source Photometry')\n",
    "plt.xlabel(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')))\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Conclusions\n",
    "\n",
    "The `photutils` package provides a comprehensive toolkit for astronomers to perform aperture photometry, including customizable aperture shapes that allow for more precise photometry and easy photometric correction.\n",
    "\n",
    "**To continue with this `photutils` tutorial, go on to the [PSF photometry notebook](../04_psf_photometry/04_psf_photometry.ipynb).**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Additional Resources\n",
    "For more examples and details, please visit the [photutils](http://photutils.readthedocs.io/en/stable/index.html) documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## About this Notebook\n",
    "**Authors:** Lauren Chambers (lchambers@stsci.edu), Erik Tollerud (etollerud@stsci.edu), Tom Wilson (towilson@stsci.edu) Clare Shanahan (cshanahan@stsci.edu)\n",
    "<br>**Updated:** May 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Page](#title_ID)\n",
    "<img style=\"float: right;\" src=\"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"STScI logo\" width=\"200px\"/>"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}