04_psf_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",
    "# PSF Photometry with `photutils`\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### What is PSF photometry?\n",
    "A more specific form of photometry than aperture photometry, PSF photometry takes into account the shape of a source's point spread function (PSF). The PSF is a model that represents the distribution of light from a point source as it falls onto a detector. An example of a basic PSF is simply a 2-D Gaussian, while more complex PSFs can include distortion, diffraction, or interference effects associated with a particular telescope. For instance, the PSFs from the Hubble Space Telescope and the James Webb Space Telescope have been meticulously modeled, and can be simulated with the [Tiny Tim](http://www.stsci.edu/hst/observatory/focus/TinyTim) and [WebbPSF](https://github.com/mperrin/webbpsf) software packages, respectively. However, for datasets that do not have readily available PSF models, such models can be statistically generated by analyzing the image itself.\n",
    "\n",
    "The `photutils` package provides tools that combine background estimation, source detection, and model-fitting to perform PSF photometry on image data.\n",
    "\n",
    "##### What does this tutorial include?\n",
    "This tutorial covers how to perform PSF photometry with `photutils`, including the following methods:\n",
    "* Gaussian PSF Photometry\n",
    "* Iterative Subtraction\n",
    "* Point Response Function (PRF) Photometry\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.psf` [documentation](http://photutils.readthedocs.io/en/stable/psf.html).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    \n",
    "**Warning:** The PSF photometry API is currently considered experimental and may change in the future. The photutils development team will aim to keep compatibility where practical, but will not finalize the API until sufficient user feedback has been accumulated.\n",
    "\n",
    "</div>\n",
    "\n",
    "<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\n",
    "from astropy.visualization import ImageNormalize, LogStretch\n",
    "import matplotlib\n",
    "from matplotlib.colors import LogNorm\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import LogLocator\n",
    "import numpy as np\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": [
    "## Circular Apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating Apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performing Aperture Photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating Aperture Corrections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elliptical Apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating Apertures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performing Aperture Photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating Aperture Corrections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Exercises"
   ]
  },
  {
   "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)\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
}