LocalCoordAperturePhotometry.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local coordinate aperture photometry for input to DSLR data worksheet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from random import random\n",
    "from collections import OrderedDict\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from astropy.io import fits\n",
    "from astropy.visualization import astropy_mpl_style\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "from matplotlib.patches import Circle\n",
    "from matplotlib.offsetbox import (TextArea, DrawingArea, OffsetImage, AnnotationBbox)\n",
    "plt.style.use(astropy_mpl_style)\n",
    "%matplotlib inline\n",
    "\n",
    "from PythonPhot import aper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Photometry of a list of FITS files, creating a table of times and instrumental magnitudes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def multi_file_photometry(fits_root, fits_files, data_index, coords, dataframe, \n",
    "                          aperture_radius, inner_sky_radius, outer_sky_radius, \n",
    "                          gain=1, zeropoint=0, suffix='.fit'):\n",
    "    \n",
    "    for fits_file in fits_files:\n",
    "        fits_file_path = os.path.join(fits_root, fits_file)\n",
    "        hdus = fits.open(fits_file_path)\n",
    "            \n",
    "        instr_mags = []\n",
    "        for x, y in coords:\n",
    "            time, mag = aperture_photometry(hdus[data_index], x, y, \n",
    "                                            aperture_radius, inner_sky_radius, outer_sky_radius,\n",
    "                                            gain=gain, zeropoint=zeropoint)\n",
    "            instr_mags.append(mag)\n",
    "        \n",
    "        dataframe[fits_file[0:fits_file.rindex(suffix)]] = [time] + instr_mags  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Single image+coordinate photometry, returning a time and instrumental magnitude. Invoked by multi_file_photometry()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def aperture_photometry(hdu, x, y, \n",
    "                        aperture_radius, inner_sky_radius, outer_sky_radius, \n",
    "                        gain=1, zeropoint=0):\n",
    "\n",
    "    image_data = hdu.data\n",
    "    time = hdu.header[time_name]\n",
    "        \n",
    "    mag, magerr, flux, fluxerr, sky, skyerr, badflag, outstr = \\\n",
    "                aper.aper(image_data, x, y, phpadu=gain, \n",
    "                          apr=aperture_radius, zeropoint=zeropoint,\n",
    "                          skyrad=[inner_sky_radius, outer_sky_radius], \n",
    "                          exact=True)\n",
    "\n",
    "    return time, mag[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Display an image with target and reference stars annotated, to sanity check local coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def show_image(image_data, coord_map, aperture_size, annotate=True, vmin=10, vmax=200, figx=20, figy=10):\n",
    "    fig = plt.figure(figsize=(figx, figy))\n",
    "    plt.imshow(image_data, cmap='gray', vmin=vmin, vmax=vmax)\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.colorbar()\n",
    "\n",
    "    if annotate:\n",
    "        for designation in coord_map:\n",
    "            xy = coord_map[designation]\n",
    "            annotate_image(fig.axes[0], designation, xy, aperture_size)\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Annotate plot axis with coordinate positions and designations. Invoked by show_image()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def annotate_image(axis, designation, xy, aperture_size):\n",
    "    axis.plot(xy[0], xy[1], 'o', markersize=aperture_size, \n",
    "              markeredgecolor='r', markerfacecolor='none', \n",
    "              markeredgewidth=2)\n",
    "    \n",
    "    offsetbox = TextArea(designation, minimumdescent=False)\n",
    "\n",
    "    ab = AnnotationBbox(offsetbox, xy,\n",
    "                        xybox=(-20, 40+random()*10-10),\n",
    "                        xycoords='data',\n",
    "                        boxcoords=\"offset points\",\n",
    "                        arrowprops=dict(arrowstyle=\"->\"))\n",
    "    \n",
    "    axis.add_artist(ab)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inputs\n",
    "### Change these to suit your environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### File settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Instrumental magnitude output file path\n",
    "instr_mag_file_root = \"/Users/david/_photometry_working\"\n",
    "instr_mag_csv_file = \"instr_mags.csv\"\n",
    "\n",
    "# FITS file directory\n",
    "fits_root = \"/Users/david/_photometry_working\"\n",
    "\n",
    "# B, G, and R FITS file prefixes to identify files,\n",
    "# e.g. stk-median-g matches stk-median-g1.fit, stk-median-g2.fit, ... \n",
    "fits_prefixes = [\"stk-median-b\", \"stk-median-g\", \"stk-median-r\"]\n",
    "\n",
    "# FITS file data HDU index\n",
    "data_index = 0\n",
    "\n",
    "# Time column name\n",
    "time_name = \"JD\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Map of object designations to local coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Ordered dictionary of object names/IDs to local coordinates\n",
    "position_map = OrderedDict({\n",
    " # ** Your name-tuple pairs go here **\n",
    " \"eta Car\":(2374.478,814.152),\n",
    " \"45 (CHK)\":(2352.016,1629.581),\n",
    " \"46\":(2585.264,600.724),\n",
    " \"47\":(2512.030,1096.505),\n",
    " \"51\":(2664.810,1040.385),\n",
    " \"52 795\":(1657.205,1312.103),\n",
    " \"52 708\":(1777.295,1584.834),\n",
    " \"55\":(2482.270,1252.879)\n",
    " # ** END **\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Aperture radii and gain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Aperture radii\n",
    "measurement_aperture = 13\n",
    "inner_sky_annulus = 16\n",
    "outer_sky_annulus = 21\n",
    "\n",
    "# ph/ADU\n",
    "# Note: PythonPhot's aperture photometry function takes a phadu parameter.\n",
    "# Assumption: this is photons/ADU or e-/ADU, i.e. gain.\n",
    "gain=1.67"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find B, G, R files in the FITS file directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "files = os.listdir(fits_root)\n",
    "\n",
    "fits_files = []\n",
    "for fits_prefix in fits_prefixes:\n",
    "    fits_files += sorted([file for file in files if fits_prefix in file])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aperture location sanity check by visual inspection\n",
    "#### Arbitrarily choose the first G FITS file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fits_file = fits_files[0]\n",
    "\n",
    "hdus = fits.open(os.path.join(fits_root, fits_file))\n",
    "image_data = hdus[data_index].data\n",
    "\n",
    "median = np.median(image_data)\n",
    "show_image(image_data, position_map, measurement_aperture, annotate=True, vmin=10, vmax=median*4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aperture photometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create empty table with time and object headers\n",
    "pd.options.display.float_format = '{:,.6f}'.format\n",
    "instr_mag_df = pd.DataFrame()\n",
    "names = [name for name in position_map]\n",
    "instr_mag_df['name'] = [time_name] + names\n",
    "instr_mag_df.set_index('name', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Carry out photometry on B, G, R FITS files, yielding instrumental magnitudes\n",
    "positions = position_map.values()\n",
    "\n",
    "multi_file_photometry(fits_root, fits_files, data_index, positions, instr_mag_df, \n",
    "                      measurement_aperture, inner_sky_annulus, outer_sky_annulus, gain=gain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save photometry table as CSV\n",
    "instr_mag_csv_path = os.path.join(instr_mag_file_root, instr_mag_csv_file)\n",
    "instr_mag_df.T.to_csv(instr_mag_csv_path)\n",
    "\n",
    "# Display photometry table\n",
    "instr_mag_df.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "* http://docs.astropy.org/en/stable/io/fits\n",
    "* http://docs.astropy.org/en/stable/generated/examples/io/plot_fits-image.html\n",
    "* http://www.astropy.org/astropy-tutorials/FITS-images.html\n",
    "* https://github.com/djones1040/PythonPhot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}