08-WCS/WCS_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WCS solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Refine a WCS using a list of detections and a reference catalog\n",
    "\n",
    "Refine a WCS for a science image exposure from the Zwicky Transient Facility from these ingredients:\n",
    "* An initial header\n",
    "* A detection list cut at 17th magnitude, in file `data/ztf_detections_17thmag.csv`\n",
    "* A reference catalog with coordinates and magnitudes from Gaia cut at 17 Gaia G magnitude, in `data/Gaia-gaia_dr2_source-ztf-20190606224213_000667_zr.csv`\n",
    "\n",
    "The exercise makes use of `astropy.wcs`, `astropy.coordinates` and the projection capabilities of WCSAxes.\n",
    "\n",
    "1. Read in the detection list and the reference catalog with `astropy.table.Table.read`\n",
    "2. Calculate starting RAs and Decs for the detection list using the initial WCS\n",
    "3. Create SkyCoord instances for the initial detection coordinates and the Gaia coordinates\n",
    "4. Plot the detection list and the Gaia list in a scatter plot\n",
    "5. Match the detection list and the Gaia list\n",
    "6. Refine the WCS using the `fit_wcs_from_points` function from `astropy.wcs.utils`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import everything we'll need for the exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Table\n",
    "from astropy.wcs import WCS\n",
    "from astropy.wcs.utils import fit_wcs_from_points\n",
    "import astropy.units as u\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the initial WCS programatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_wcs = WCS(naxis=2)\n",
    "initial_wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']  \n",
    "initial_wcs.wcs.crval = [149.07662386535503, 33.32164150821777]  \n",
    "initial_wcs.wcs.crpix = [-3305.678, -7136.481]\n",
    "initial_wcs.wcs.cd = [[-0.0002817188, -1.554e-07],\n",
    "                      [-1.998e-07, -0.0002819204]]  \n",
    "initial_wcs.array_shape = [3080, 3072] # NAXIS2, NAXIS1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_wcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Read in the detection list and the reference catalog\n",
    "\n",
    "Read in the detections and the reference catalog using `astropy.table.Table` with `format='csv'`.\n",
    "The detections table is in `'data/ztf_detections_17thmag.csv'` and the reference catalog is `'data/Gaia-gaia_dr2_source-ztf-20190606224213_000667_zr.csv'`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detections = Table.read(os.path.join('data', 'ztf_detections_17thmag.csv'), format='csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_catalog = Table.read(os.path.join('data', 'Gaia-gaia_dr2_source-ztf-20190606224213_000667_zr.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Calculate starting RAs and Decs for the detection list using the initial WCS\n",
    "\n",
    "Use the `initial_wcs.all_pix2world` function to calculate starting RA and Dec from the `detections['xpos']` and `detections['ypos']` columns. The pixel positions use the FITS numbering convention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_ra, initial_dec = initial_wcs.all_pix2world(detections['xpos'], detections['ypos'], 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Create SkyCoord instances for the initial detection coordinates and the Gaia coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_coords = SkyCoord(ra=initial_ra, dec=initial_dec, unit=u.deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_coords = SkyCoord(ra=ref_catalog['ra'], dec=ref_catalog['dec'], unit=u.deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Plot the detection list and the Gaia list in a scatter plot\n",
    "\n",
    "Use `projection=initial_wcs` to make a scatter plot using `gaia_coords` and `initial_coords`. The open circles are sized according to magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,8))\n",
    "ax = plt.subplot(projection=initial_wcs)\n",
    "\n",
    "ax.scatter(gaia_coords.ra, \n",
    "           gaia_coords.dec,  c=None, marker='o',\n",
    "           s=20*(18 - ref_catalog['phot_g_mean_mag']),\n",
    "           facecolors='none', edgecolors='green',\n",
    "           transform=ax.get_transform('world'))\n",
    "\n",
    "ax.scatter(initial_coords.ra, \n",
    "           initial_coords.dec, c=None, marker='o',\n",
    "           s=20*(18 - detections['mag']),\n",
    "           facecolors='none', edgecolors='blue',\n",
    "           transform=ax.get_transform('world'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Match the detection list and the Gaia list\n",
    "\n",
    "Use the `initial_coords.search_around_sky` method with a 15 arcsecond radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idxgaia, idxdet, d2d, d3d = initial_coords.search_around_sky(gaia_coords, 15*u.arcsec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_matched = gaia_coords[idxgaia]\n",
    "detections_xpos_matched = detections['xpos'][idxdet]\n",
    "detections_ypos_matched = detections['ypos'][idxdet]\n",
    "print(len(gaia_matched), len(detections_xpos_matched), len(detections_ypos_matched))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Refine the WCS using the `fit_wcs_from_points` function\n",
    "\n",
    "Look at the help for `fit_wcs_from_points` and use it to fit a new WCS.\n",
    "\n",
    "Optionally, calculate new RAs and Decs for the matched pixel coordinates, and make another scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_wcs_from_points?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_wcs = fit_wcs_from_points((detections_xpos_matched, detections_ypos_matched),\n",
    "                                 gaia_matched,\n",
    "                                 projection='TAN', sip_degree=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_wcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the SIP distortion coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_wcs.sip.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_wcs.sip.b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optionally, calculate new RAs and Decs for the matched pixel coordinates, and make another scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_ra, fitted_dec = fitted_wcs.all_pix2world(detections_xpos_matched,\n",
    "                                                 detections_ypos_matched, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,8))\n",
    "ax = plt.subplot(projection=fitted_wcs)\n",
    "\n",
    "ax.scatter(gaia_matched.ra, \n",
    "           gaia_matched.dec,  c=None, marker='o',\n",
    "           s=20*(18 - ref_catalog['phot_g_mean_mag'][idxgaia]),\n",
    "           facecolors='none', edgecolors='green',\n",
    "           transform=ax.get_transform('world'))\n",
    "\n",
    "ax.scatter(fitted_ra, \n",
    "           fitted_dec,  c=None, marker='o',\n",
    "           s=20*(18 - detections['mag'][idxdet]),\n",
    "           facecolors='none', edgecolors='blue',\n",
    "           transform=ax.get_transform('world'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}