09-Photutils/03-aperture_local_bkgsub.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"data/photutils_banner.svg\" width=500 alt=\"Photutils logo\" style=\"margin-left: 0;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<h2 style=\"margin-top: 0\">In this notebook, we will cover:</h2>\n",
    "\n",
    "-  Aperture photometry with local-background subtraction using CircularAnnulus apertures\n",
    "  - Using the mean within a circular annulus\n",
    "  - Using the sigma-clipped median within a circular annulus\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local Background Subtraction"
   ]
  },
  {
   "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/**](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",
    "from astropy.wcs import WCS\n",
    "\n",
    "sci_fn = 'data/xdf_hst_wfc3ir_60mas_f160w_sci.fits'\n",
    "rms_fn = 'data/xdf_hst_wfc3ir_60mas_f160w_rms.fits'\n",
    "\n",
    "hdr = fits.getheader(sci_fn)\n",
    "data = fits.getdata(sci_fn)\n",
    "error = fits.getdata(rms_fn)\n",
    "wcs = WCS(hdr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the total error:  background error plus source Poisson error\n",
    "from photutils.utils import calc_total_error\n",
    "\n",
    "eff_gain = hdr['TEXPTIME']\n",
    "total_error = calc_total_error(data, error, eff_gain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The background in the XDF image has already been subtracted.  Let's add a background of `7.2 e-/s` (per pixel) so we'll have something to subtract in these examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data += 7.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# display the data\n",
    "from astropy.visualization import simple_norm\n",
    "\n",
    "plt.figure(figsize=(5, 5))\n",
    "norm = simple_norm(data, 'sqrt', percent=99.5)\n",
    "plt.imshow(data, norm=norm)\n",
    "cbar = plt.colorbar(shrink=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Perform aperture photometry at multiple positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from photutils.aperture import CircularAperture\n",
    "\n",
    "# define three apertures\n",
    "positions = [(90.73, 59.43), (73.63, 139.41), (43.62, 61.63)]\n",
    "radius = 5.0  # pixels\n",
    "circ_aper = CircularAperture(positions, r=radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the apertures\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.imshow(data, norm=norm)\n",
    "circ_aper.plot(color='red', lw=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# perform aperture photometry for the three sources\n",
    "# data here includes background, so the aperture sums are *not* the source fluxes\n",
    "# ideally the data array should be background subtracted before performing aperture photometry\n",
    "import astropy.units as u\n",
    "from photutils.aperture import ApertureStats\n",
    "\n",
    "unit = u.electron / u.s\n",
    "circ_apstats = ApertureStats(data << unit, circ_aper, error=total_error << unit)\n",
    "circ_apstats.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local background estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One often wants to also estimate the local background around each source using a nearby aperture or annulus aperture surrounding each source.  A simple method for doing this is to use the `ApertureStats` class to compute the mean background level within the background aperture. This class can also be used to calculate more advanced statistics (e.g., a sigma-clipped median) within the background aperture (e.g., a circular annulus). We show examples of both below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by creating circular-annulus apertures at the same positions as the `CircularApertures` defined above. We'll define each circular annulus to have an inner and outer radius of 10 and 15 pixels, respectively.\n",
    "\n",
    "The circular-annulus apertures will be used for local background estimation around the sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils.aperture import CircularAnnulus\n",
    "\n",
    "bkg_aper = CircularAnnulus(positions, r_in=10.0, r_out=15.0)\n",
    "apers = [circ_aper, bkg_aper]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the apertures\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.imshow(data, norm=norm)\n",
    "aper_patches = circ_aper.plot(color='white', lw=2, label='Photometry aperture')\n",
    "bkg_patches = bkg_aper.plot(color='orange', lw=2, label='Background annulus', hatch='///')\n",
    "handles = (aper_patches[0], bkg_patches[0])\n",
    "plt.legend(loc=(0.03, 0.85), facecolor='#458989', labelcolor='white',\n",
    "           handles=handles, prop={'weight': 'bold', 'size': 11});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple mean within a circular annulus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This first example uses the mean value in circular annulus as the background value.\n",
    "\n",
    "We'll use the `ApertureStats` class to compute the mean background level (per pixel) within the annulus aperture at each position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from photutils.aperture import ApertureStats\n",
    "\n",
    "bkg_apstats = ApertureStats(data << unit, bkg_aper)\n",
    "bkg_mean = bkg_apstats.mean\n",
    "print(bkg_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total background within the circular apertures is then the per-pixel background level multiplied by the circular-aperture areas.  The circular-aperture areas can be calculated from the `ApertureStats` object with a circular aperture (defined above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_bkg = bkg_mean * circ_apstats.sum_aper_area.value\n",
    "print(total_bkg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's calculate the local background-subtracted sum within the circular apertures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apersum_bkgsub = circ_apstats.sum - total_bkg\n",
    "print(apersum_bkgsub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sigma-clipped median within a circular annulus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, the local background level around each source is estimated as the sigma-clipped median value within the circular annulus. We’ll use the `ApertureStats` class to compute both the photometry (aperture sum) and the background level.\n",
    "\n",
    "We'll use the same circular and circular-annulus apertures defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.stats import SigmaClip\n",
    "\n",
    "sigclip = SigmaClip(sigma=3.0, maxiters=10)\n",
    "circ_stats = ApertureStats(data << unit, circ_aper, sigma_clip=None)\n",
    "bkg_stats = ApertureStats(data << unit, bkg_aper, sigma_clip=sigclip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's calculate the sigma-clipped median values in the background annulus apertures.\n",
    "\n",
    "Note that the values are *very* close to the background value of `7.2` that we added to the data at the beginning of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_stats.median"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, the total background within the circular apertures is then the per-pixel background level multiplied by the circular-aperture areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_bkg = bkg_stats.median * circ_stats.sum_aper_area.value\n",
    "print(total_bkg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's calculate the local background-subtracted sum within the circular apertures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "apersum_bkgsub = circ_stats.sum - total_bkg\n",
    "print(apersum_bkgsub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if you want to compute all of the source properties (i.e., in addition to only sum) on the local-background-subtracted data, you may input the per-pixel local background values to `ApertureStats` via the `local_bkg` keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aper_stats_bkgsub = ApertureStats(data << unit, circ_aper,\n",
    "                                  local_bkg=bkg_stats.median << unit)\n",
    "print(aper_stats_bkgsub.sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note these background-subtracted values are the same as those above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning alert-block\">\n",
    "<h3 style='margin-top: 0;'>Learn More</h3>\n",
    "\n",
    "The [segmentation notebook](04-segmentation.ipynb) covers:\n",
    "    \n",
    "- The basics of image segmentation\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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}