PSF_pyinthesky_API.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy.nddata import NDData\n",
    "from photutils import psf, background, daofind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from stsci.jwst.io import load_nircam\n",
    "from stsci.jwst import jwst_phot_tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load up the data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "program = load_nircam('my_cycle1_data_because_all_my_friends_were_on_the_TAC.asdf')\n",
    "images = program[program.filter == 'F070W'].images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`program` is something that behaves like an association table, but has a fairly easy-to-use interface for selecting individual exposures, like how it's used here.\n",
    "\n",
    "`images` is a *list* of [NDData](http://docs.astropy.org/en/stable/api/astropy.nddata.NDData.html) objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stacked_image = program.drizzled_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the above I'm assuming there's a pipeline step that does combining on stacks in a reasonably sensible way (which I'm speculating to be drizzle-like.  Many science users will want to do that themselves, which of course they should be able to do if they want to (presumably inside this notebook).  But either way, `stacked_image` should come out as a single `NDData` object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inital prep work (as needed) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "full_bkgs = [background.Background(image, ...) for image in images]\n",
    "images_sub = [image - full_bkg for img, full_bkg in zip(images, full_bkgs)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stacked_bkg = background.Background(stacked_image, ...)\n",
    "stacked_sub = stacked_image - stacked_bkg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It might need to be a *bit* more complicated than the above to get science-quality subtraction... but ideally something nearly that simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the output of `load_nircam` might yield an *already-subtracted* `NDData`... and also a `background.Background` object that the user can add back in of they want the non-subtracted version."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the PSF, or derive it if needed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The \"easy\" way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psf_model = jwst_phot_tools.get_psf(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`psf_model` above is some kind of astropy.modeling.models 2D model.  It follows the interface described at https://photutils.readthedocs.io/en/latest/api/photutils.psf.psf_photometry.html#photutils.psf.psf_photometry\n",
    " \n",
    "`image` already contains the metadata (equivalent to the fits header keywords), and `get_psf` does some kind of fancy black magic based on that + an optical model of JWST or whatever to figure out what the PSF is. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The \"hard\" way "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method identifies and then fits PSF stars.  It does the finding on the stacked image, but fitting PSFs on individual exposures.  Detection could also be done on individual exposures, but that just adds a bunch of for loops here, which is a bit harder to read."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from photutils import daofind\n",
    "from photutils import aperture_photometry\n",
    "from astropy.stats import sigma_clipped_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find sources to be possible PSF stars - for this example we'll use daofind, but whatever could be given "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mean, median, std = sigma_clipped_stats(combined_image.data, sigma=3.0, iters=5) \n",
    "sources = daofind(combined_image, fwhm=5, threshold=5.*std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A slightly more convenient version of the above cell, not yet supported in the current photutils:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sources = daofind(stacked_sub, fwhm=0.1*u.arcsec, threshold=5.*std)  \n",
    "# the fwhm argument would know to use `stacked_sub.wcs` to determine the pixel scale - daofind doesn't do that yet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Do quick aperture photometry "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "positions = [(row['xcentroid'], row['ycentroid']) for row in sources]\n",
    "apertures = CircularAperture(positions, r=15) \n",
    "ap_phot = aperture_photometry(stacked_sub, apertures)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also is not-yet-supported (?), but probably should be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "apertures = CircularAperture(sources, r=15)\n",
    "# CircularAperture currently doesn't know how to swallow tables that come out of the star-finders, but should\n",
    "ap_phot = aperture_photometry(stacked_sub, apertures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# maybe this should autmoatically go into aperture_photometry?\n",
    "ap_phot['instrumental_mags'] = u.Magnitude(ap_phot['aperture_sum']*u.count)\n",
    "# or maybe JWST will provide:\n",
    "ap_phot['calibrated_mags'] = jwst_phot_tools.calibrate(ap_phot['aperture_sum'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now pick out the PSF stars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psf_stars = psf.pick_stars(ap_phot, image_sub, nstars=10, min_seperation=50*u.pixel, psfrng=(10*u.mag, 12*u.mag))\n",
    "# could give `psfrng` in units of counts, if desired"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`psf_stars` is a Table that is a subset of `ap_phot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for star in psf_stars:\n",
    "    # this already exists in NDData\n",
    "    cutout = nddata.Cutout2D(image_sub.data, (star['xcentroid'], star['ycentroid']), size=100*u.pixel)\n",
    "    plt.imshow(cutout.data)\n",
    "    \n",
    "# might want to provide a simple quick function to do the same as the above along the lines of:\n",
    "psf.show_psf_cutouts(psf_stars, image_sub, size=100*u.pixel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't like any of them, just eliminate them from the table like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "del psf_stars[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### And build the PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psf_model = jwst_phot_tools.NIRCAMPSF.make(psf_stars, images_sub)\n",
    "# this would yield whatever kind of PSF the NIRCAM folks decide is \"best\" - perhaps a subclass of Jay's PSF model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or a user might prefer their own model for the PSF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psf_model0 = picky_users_personal_tools.SomeFancy2DModel()\n",
    "psf_model = picky_users_personal_tools.psf_fitter(psf_model0, psf_stars, images_sub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The \"intermediate\" way "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psf_model = psf.create_psf(images_sub, find_fwhm=0.1*u.arcsec, find_threshold='5sigma', \n",
    "                           aperture_size=0.3*u.arcsec, psftype=jwst_phot_tools.JWSTPSF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Would do all the steps as above, but just in a single easy-to-call function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now do the actual PSF photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forced photometry "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This assumes you went the hard way, which has done source-finding already, using the coadded/stacked image. Effectively that's doing forced-photometry, which is conceptually straightforward, but requires that the catalogs be  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the \"mode\" currently only supports \"sequential\", but presumably something like NStar will work in the near future\n",
    "psf_phot_tables = []\n",
    "for image_sub in images_sub:\n",
    "    psf_phot_table = psf.psf_photometry(image_sub, sources, mode='nstar', model=psf_model)\n",
    "    psf_phot_tables.append(psf_phot_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe you also want local background subtraction? (this is not yet in the function):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bkg = background.MMMBackground(iters=5)  # for some reason you don't trust it for the full 20 iters\n",
    "bkg_region = CircularAnnulus(r_in=15, r_out=30)  \n",
    "\n",
    "psf_phot_tables = []\n",
    "for image_sub in images_sub:\n",
    "    psf_phot_table = psf.psf_photometry(image_sub, sources, mode='nstar', model=psf_model, \n",
    "                                        background_sub=bkg, background_aperture=bkg_region)\n",
    "    psf_phot_tables.append(psf_phot_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psf_phot_table = psf.combine_phot_tables(psf_phot_tables)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function above takes the single-exposure tables and adds up all the `flux_fit`'s for matching objects.  Lots of prior art in that, but also lots of ways to go wrong.  This is possibly better for some ground-based datasets (where distortion and such can change a lot from exposure-to-exposure), but is probably not good for cases where you can do simultaneous matching.  It's probably fine (better?) for reasonably high S/N stars, though."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simultaneous fitting "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the telescope is fairly stable over multiple exposures/dithers, it's probably best to instead fit multiple images *simultaneously*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "psf_phot_table = []\n",
    "psf_phot_table = psf.psf_photometry(images_sub, sources, model=psf_model, mode='nstar', )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "this does not currently work in `psf_photometry`, but it *could*  be changed to interpret \"list of NDData\" as  indicating that simultaneous fitting is requested.  The resulting table would look just like the single-fit version, but of course the actual fit process would be very different.\n",
    "\n",
    "In principal it does *not* require a separate fitter/model, however.  Instead, the fitter should be given *multiple* images and matching `x` values in one go.  This requires careful documentation of how the models are written, though: e.g., the centers need to be in \"world\" coordinates, even though the fitting is done in pixel coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next is the science "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "calibrated_mags = jwst_phot_tools.calibrate(psf_phot['flux_fit'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.hist(calibrated_mags)  # And there's your luminosity function!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "win_nobel_prize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}