09-Photutils/01-aperture_basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"data/photutils_banner.svg\" width=500 alt=\"Photutils logo\" style=\"margin-left: 0;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<h2 style=\"margin-top: 0\">In this notebook, we will cover:</h2>\n",
    "\n",
    "- The basics to perform aperture photometry with the `aperture_photometry` function\n",
    "- Aperture photometry (and other properties) with the `ApertureStats` class\n",
    "\n",
    "This notebook builds on the previous tutorials for Astropy Units/Quantities, Coordinates, FITS, and Tables.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Initial imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Change some default plotting parameters\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['image.origin'] = 'lower'\n",
    "mpl.rcParams['image.interpolation'] = 'nearest'\n",
    "\n",
    "# Run the %matplotlib magic command to enable inline plotting\n",
    "# in the current notebook.  Choose one of these:\n",
    "%matplotlib inline\n",
    "# %matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We'll start by reading science data and error arrays from FITS files located in the **data/** subdirectory.  The FITS files contain 2D cutout images from the [Hubble Extreme-Deep Field (XDF)](https://archive.stsci.edu/prepds/xdf/) taken with the [Wide Field Camera 3 (WFC3)](https://www.stsci.edu/hst/instrumentation/wfc3) IR channel in the F160W filter (centered at ~1.6 $\\mu m$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "sci_fn = 'data/xdf_hst_wfc3ir_60mas_f160w_sci.fits'\n",
    "rms_fn = 'data/xdf_hst_wfc3ir_60mas_f160w_rms.fits'\n",
    "sci_hdulist = fits.open(sci_fn)\n",
    "rms_hdulist = fits.open(rms_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print some information about the science data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sci_hdulist.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data array is in the **`0`** extension, has a shape of `200 x 200` pixels, and has a data type of `float32` (32 bit floating-point numbers).\n",
    "\n",
    "Let's extract the data and error arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "data = sci_hdulist[0].data.astype(float)\n",
    "error = rms_hdulist[0].data.astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`data` and `error` are now 2D [numpy](https://numpy.org) arrays, each with a shape of 200 x 200 pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape, error.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's extract the data header and create an Astropy [World Coordinate System (WCS)](https://docs.astropy.org/en/stable/wcs/index.html) object from the FITS header WCS information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.wcs import WCS\n",
    "\n",
    "hdr = sci_hdulist[0].header\n",
    "wcs = WCS(hdr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's display the science image.  Here we use the [astropy.visualization](https://docs.astropy.org/en/stable/visualization/index.html) subpackage to apply a square-root stretch to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.visualization import simple_norm\n",
    "\n",
    "plt.figure(figsize=(5, 5))\n",
    "norm = simple_norm(data, 'sqrt', percent=99.0)\n",
    "plt.imshow(data, norm=norm)\n",
    "plt.title('HST F160W XDF Cutout');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Aperture Objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Photutils provides circular, elliptical, and rectangular aperture shapes (plus annulus versions of each).\n",
    "\n",
    "Further, there are two types of aperture classes, defined either with pixel or sky (celestial) coordinates.\n",
    "\n",
    "These are the names of the [aperture classes](https://photutils.readthedocs.io/en/stable/user_guide/aperture.html#apertures) that are defined in pixel coordinates:\n",
    "\n",
    "* `CircularAperture`\n",
    "* `CircularAnnulus`\n",
    "\n",
    "* `EllipticalAperture`\n",
    "* `EllipticalAnnulus`\n",
    "\n",
    "* `RectangularAperture`\n",
    "* `RectangularAnnulus`\n",
    "\n",
    "The aperture classes defined in celestial coordinates have `Sky` prepended to their names:\n",
    "\n",
    "* `SkyCircularAperture`\n",
    "* `SkyCircularAnnulus`\n",
    "\n",
    "* `SkyEllipticalAperture`\n",
    "* `SkyEllipticalAnnulus`\n",
    "\n",
    "* `SkyRectangularAperture`\n",
    "* `SkyRectangularAnnulus`\n",
    "\n",
    "\n",
    "The aperture shapes look like this:\n",
    "\n",
    "<img src='data/apertures.png' alt='Figure of aperture shapes' width=700px>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Methods for handling aperture/pixel overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The apertures will usually only partially overlap some of the pixels in the data.\n",
    "\n",
    "Photutils apertures provide three methods for handling the aperture overlap with the pixel grid of the data array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"data/photutils_aperture_methods.svg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the default method (`method='exact'`), the exact fractional overlap of the aperture with each pixel is calculated. Each pixel value is then weighted by the exact fractional overlap of the aperture.\n",
    "\n",
    "For `method='center'`, a pixel is considered to be entirely in or out of the aperture depending on whether its center is located within the aperture.\n",
    "\n",
    "For `method='subpixel'`, pixels are divided into a number of subpixels, which are in or out of the aperture based on their centers. Each pixel value is then weighted by its fraction of subpixel overlaps. For this method, the number of subpixels needs to be set with the `subpixels` keyword. The `subpixels` keyword is ignored for the `'exact'` and `'center'` methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating Aperture Objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's define a circular aperture at a given (x, y) pixel position and radius (in pixels).\n",
    "\n",
    "Photutils has many tools for detecting sources (e.g., see [photutils.detection](https://photutils.readthedocs.io/en/stable/user_guide/detection.html) and [photutils.segmentation](https://photutils.readthedocs.io/en/stable/user_guide/segmentation.html)) and measuring source positions (e.g., [photutils.centroids](https://photutils.readthedocs.io/en/stable/user_guide/centroids.html)). For this example, we'll assume we've already measured the source centroid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from photutils.aperture import CircularAperture\n",
    "\n",
    "position = (90.73, 59.43)  # (x, y) centroid position (pixels)\n",
    "radius = 5.0  # pixels\n",
    "aperture = CircularAperture(position, r=radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `aperture` variable is now an instance of a [CircularAperture](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAperture.html) object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aperture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the aperture on the data using its `plot` method, with a custom color and line width."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "plt.imshow(data, norm=norm)\n",
    "aperture.plot(color='red', lw=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performing Aperture Photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's perform photometry on the XDF data using this circular aperture.\n",
    "\n",
    "After the aperture object is created, we can  perform the photometry using the [aperture_photometry](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.aperture_photometry.html) function or the [ApertureStats](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.ApertureStats.html) class.  We'll start by using the `aperture_photometry` function.  The `ApertureStats` class will be introduced later in this notebook.\n",
    "\n",
    "In the most basic case, we simply need to input the data array and the aperture object. The default aperture overlap method is `'exact'`.\n",
    "\n",
    "Note that the input data array is assumed to be background subtracted.  If that is not the case, please see the documentation for the [photutils.background](https://photutils.readthedocs.io/en/stable/user_guide/background.html) subpackage for tools to help subtract the background.\n",
    "\n",
    "The background has already been subtracted from our XDF example data, so we skip that step in this example.\n",
    "\n",
    "<div class=\"alert alert-warning alert-block\">\n",
    "<h3 style=\"margin-top:0;\">Learn More</h3>\n",
    "\n",
    "See the [local background subtraction notebook](03-aperture_local_bkgsub.ipynb) for examples of aperture photometry with local background subtraction.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from photutils.aperture import aperture_photometry\n",
    "\n",
    "phot = aperture_photometry(data, aperture)\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is an Astropy `QTable` (Quantity Table) with several columns.  The sum of the data values within the circular aperture is in the `aperture_sum` column.\n",
    "\n",
    "The table also contains metadata, which is accessed by the `meta` attribute of the table.  The metadata is stored in a python dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot.meta['aperture_photometry_args']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other aperture overlap methods are specified using the `method` keyword.  Aperture photometry using the **'center'** method gives a slightly less accurate answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "phot = aperture_photometry(data, aperture, method='center')\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now perform aperture photometry using the `'subpixel'` method with `subpixels=5`:\n",
    "\n",
    "These parameters are equivalent to SourceExtractor aperture photometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "phot = aperture_photometry(data, aperture, method='subpixel', subpixels=5)\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Photometric Errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also input an error array to get the photometric errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "phot = aperture_photometry(data, aperture, error=error)\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The error keyword expects the **total** error array (i.e., all sources of error).\n",
    "\n",
    "However, the error array in our XDF FITS file represents only the background error.  Thus, we need to also include the source Poisson error as follows:\n",
    "\n",
    "$\\sigma_{\\mathrm{tot}} = \\sqrt{\\sigma_{\\mathrm{b}}^2 + \\frac{I}{g}}$\n",
    "  \n",
    "where $\\sigma_{\\mathrm{tot}}$ is the total error (in the same units as the data), $\\sigma_{\\mathrm{b}}$ is the background-only error, $I$ are the data values, and $g$ is the \"effective gain\".\n",
    "\n",
    "The \"effective gain\" is the value (or an array if it's variable across an image) needed to convert the data image to count units (e.g., electrons or photons), where Poisson statistics apply.\n",
    "\n",
    "Photutils provides a [calc_total_error](https://photutils.readthedocs.io/en/stable/api/photutils.utils.calc_total_error.html) function to perform the calculation of combining background-only error with the source Poisson error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this time include the Poisson error of the source\n",
    "from photutils.utils import calc_total_error\n",
    "\n",
    "# our data array is in units of e-/s\n",
    "# so the \"effective gain\" should be the exposure time\n",
    "eff_gain = hdr['TEXPTIME']  # exposure time from the FITS header\n",
    "total_error = calc_total_error(data, error, eff_gain)\n",
    "\n",
    "phot = aperture_photometry(data, aperture, error=total_error)\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aperture photometry error increased only slightly because this is a small faint source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Quantity` inputs for the `data` and `error` arrays are also allowed.  Note that the unit must be the same for the `data` and `error` inputs.\n",
    "\n",
    "The XDF dataset has units of electron/s. One way to check this is to look at the value of the `BUNIT` FITS header keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdr['BUNIT']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "\n",
    "unit = u.electron / u.s  # unit for the data and error arrays\n",
    "phot = aperture_photometry(data << unit, aperture, error=total_error << unit)\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table columns now contain `Quantity` arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot['aperture_sum']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Performing aperture photometry for multiple sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's perform aperture photometry for three sources (all with the same aperture size).  We simply define three (x, y) positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "positions = [(90.73, 59.43), (73.63, 139.41), (43.62, 61.63)]\n",
    "radius = 5.0  # pixels\n",
    "apertures = CircularAperture(positions, r=radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that an aperture object with multiple positions can be indexed or sliced to get a subset of apertures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apertures[1]  # aperture for the second position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apertures[0:2]  # apertures for the first two positions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot all three apertures on the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "plt.imshow(data, norm=norm)\n",
    "apertures.plot(color='red', lw=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's perform aperture photometry for these three sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "phot = aperture_photometry(data << unit, apertures, error=total_error << unit)\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each source is a row in the table and is given a unique **`id`** number (the first column)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we input arrays with units, the output columns are `Quantity` arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot['aperture_sum']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Including Sky Coordinates in the output table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output table will include the sky coordinates of the input x/y positions if we pass in the `WCS` transform defined when we loaded the data. Note that this requires photutils version >= 2.1.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot = aperture_photometry(data << unit, apertures, error=total_error << unit, wcs=wcs)\n",
    "phot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot['sky_center']  # a SkyCoord object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Adding columns to the photometry table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the signal-to-noise (SNR) ratio of our sources and add it as a new column to the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "snr = phot['aperture_sum'] / phot['aperture_sum_err']  # units will cancel\n",
    "phot['snr'] = snr\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above, this dataset was taken with the HST WFC3/IR F160W filter.  Using the [F160W zero point](https://spacetelescope.github.io/hst_notebooks/notebooks/WFC3/zeropoints/zeropoints.html), let's calculate the F160W AB magnitude and add it as a new column in the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "f160w_zpt = 25.941  # HST/WFC3 F160W ABmag zero point for flux in e-/s\n",
    "\n",
    "# NOTE that the log10 function can be applied only to dimensionless quantities,\n",
    "# so we use the value attribute for Quantity objects to remove the units of the aperture sum\n",
    "abmag = -2.5 * np.log10(phot['aperture_sum'].value) + f160w_zpt\n",
    "\n",
    "phot['f160w_abmag'] = abmag\n",
    "phot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving a photometry table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we write the table to an ASCII file using the ECSV format, we can read it back in later preserving all of the table metadata and `Quantity` and `SkyCoord` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aper_filename = 'f160w_aperture_photometry.ecsv'\n",
    "phot.write(aper_filename, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's read the table into a new variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import QTable\n",
    "\n",
    "tbl = QTable.read(aper_filename)\n",
    "tbl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the table metadata was preserved\n",
    "tbl.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# as was the Quantity arrays, e.g.:\n",
    "tbl['aperture_sum']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and the SkyCoord array\n",
    "tbl['sky_center']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Aperture photometry using Sky apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aperture photometry can also be performed using `SkyCoord` positions instead of pixel positions.\n",
    "\n",
    "First, let's define the sky coordinates by converting our pixel coordinates using our WCS transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "positions = [(90.73, 59.43), (73.63, 139.41), (43.62, 61.63)]\n",
    "x, y = np.transpose(positions)\n",
    "coord = wcs.pixel_to_world(x, y)\n",
    "coord"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define circular apertures centered at these sky coordinates.\n",
    "\n",
    "For sky apertures, the aperture radius must be a `Quantity` in angular units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from photutils.aperture import SkyCircularAperture\n",
    "\n",
    "radius = 0.5 * u.arcsec  # angular units\n",
    "sky_apers = SkyCircularAperture(coord, r=radius)\n",
    "sky_apers.r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When performing photometry using a sky aperture, `aperture_photometry` needs the WCS transformation, which is input via the `wcs` keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "phot_sky = aperture_photometry(data, sky_apers, wcs=wcs)\n",
    "phot_sky"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ApertureStats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [ApertureStats](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.ApertureStats.html) class can be used to create a catalog of statistics and properties for pixels within an aperture, including aperture photometry. \n",
    "\n",
    "It can calculate many properties, including statistics like `min`, `max`, `mean`, `median`, `std`, `sum_aper_area`, and `sum`. It also can be used to calculate morphological properties like `centroid`, `fwhm`, `semimajor_sigma`, `semiminor_sigma`, `orientation`, and `eccentricity`. Please see `ApertureStats` for the the complete list of properties that can be calculated. The properties can be accessed using `ApertureStats` attributes or output to an Astropy `QTable` using the `to_table` method.\n",
    "\n",
    "Most of the source properties are calculated using the \"center\" aperture-mask method, which gives aperture weights of 0 or 1. This avoids the need to compute weighted statistics for partial pixel overlaps -- the data pixel values are directly used.\n",
    "\n",
    "The `sum_method` and `subpixels` keywords are used to determine the aperture-mask method when calculating the sum-related properties: `sum`, `sum_error`, `sum_aper_area`, `data_sumcutout`, and `error_sumcutout`. The default is `sum_method='exact'`, which produces exact aperture-weighted photometry.\n",
    "\n",
    "The optional `local_bkg` keyword can be used to input the per-pixel local background of each source, which will be subtracted before computing the aperture statistics.\n",
    "\n",
    "The optional `sigma_clip` keyword can be used to sigma clip the pixel values before computing the source properties. This keyword could be used, for example, to compute a sigma-clipped median of pixels in an annulus aperture to estimate the local background level.\n",
    "\n",
    "Note that like `aperture_photometry`, `ApertureStats` expects the input data to be background subtracted.  The background has already been subtracted from our XDF example data, so we skip that step in this example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we start by creating an `Aperture` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "positions = [(90.73, 59.43), (73.63, 139.41), (43.62, 61.63)]\n",
    "radius = 5.0  # pixels\n",
    "apertures = CircularAperture(positions, r=radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create an `ApertureStats` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils.aperture import ApertureStats\n",
    "\n",
    "apstats = ApertureStats(data << unit, apertures, error=total_error << unit)\n",
    "apstats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The source properties can be accessed as attributes, for example the source sum (and errors) for aperture photometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apstats.sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apstats.sum_err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these values are exactly the same as those calculated earlier using `aperture_photometry`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phot['aperture_sum'], phot['aperture_sum_err']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Going further, here are the calculated source centroids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apstats.xcentroid, apstats.ycentroid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the source properties can be output as an astropy `QTable` using the `to_table` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apstats.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measure the circular-aperture fluxes (and errors) of the three sources centered at (x, y) positions of `(34.61, 137.99), (184.53, 157.52), and (88.75, 184.22)` within a radius of `4.5` pixels.\n",
    "\n",
    "Calculate their signal-to-noise ratios and F160W AB magnitudes.\n",
    "\n",
    "Plot the apertures on the data.\n",
    "\n",
    "Hints:\n",
    "\n",
    "* The HST/WFC3 F160W ABmag zero point for data in units of e-/s is `25.941`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you prefer to load the solution, for each cell below uncomment the cell and **run it twice** (once to load the solution and once more to run it)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load 01-aperture_exercise_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load 01-aperture_exercise_solution_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning alert-block\">\n",
    "<h3 style='margin-top: 0;'>Learn More</h3>\n",
    "\n",
    "The [aperture photometry extended notebook](02-aperture_extended.ipynb) covers:\n",
    "\n",
    "- Aperture Photometry\n",
    "    - Performing aperture photometry at multiple positions with the same aperture\n",
    "      - Using the aperture_photometry function\n",
    "      - Using the ApertureStats class\n",
    "    - Performing aperture photometry at multiple positions using different apertures\n",
    "      - Using the aperture_photometry function\n",
    "      - Using the ApertureStats class\n",
    "    - Bad pixel masking\n",
    "      - Using the aperture_photometry function\n",
    "      - Using the ApertureStats class\n",
    "    - Encircled flux\n",
    "      - Using the aperture_photometry function\n",
    "      - Using the ApertureStats class\n",
    "\n",
    "- ApertureMask objects\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}