08-WCS/WCS_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## astropy.wcs\n",
    "Implements the FITS WCS standard and some commonly used distortion conventions.\n",
    "\n",
    "This tutorial will show how to create a WCS object from a FITS file and how to use it to transform coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.utils.data import get_pkg_data_filename\n",
    "from astropy.io import fits\n",
    "from astropy import wcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a WCS object from the header of a FITS file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open a file with `astropy.fits` and look at it. (This file contains only the header and no data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sip_file_name = get_pkg_data_filename('data/sip.fits', package='astropy.wcs.tests')\n",
    "\n",
    "sip_file = fits.open(sip_file_name)\n",
    "sip_file.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a WCS object pass the header with the WCS keywords to astropy.wcs.WCS. In this case it is the primary header. A warning is issued due to the lack of image data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = wcs.WCS(sip_file[0].header)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transforming between pixel coordinates and sky coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform the transformation from detector to sky, including distortions, pass x and y and an 'origin'. The third argument, 'origin', indicates whether the coordinates are 1-based (like FITS), or 0-based (like python).\n",
    "\n",
    "The inputs can be numbers, numpy arrays or array like objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra, dec = w.all_pix2world([1, 100], [2, 200], 1)\n",
    "print(ra, dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perfom the inverse transformation - from sky to detector coordinates.\n",
    "\n",
    "If analytical inverse is not available (often the case in the presence of distortion), then an iterative inverse is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(w.all_world2pix(ra, dec, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases it is useful to omit the distortion and perform the core WCS transforms only:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra, dec = w.wcs_pix2world([1, 100], [2, 200], 1)\n",
    "print(ra, dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.wcs_world2pix(ra, dec, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a WCS programatically"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A WCS object can be created programatically. Here is a concise example with 1 arcsecond pixels that is aligned with \"North up, East to the left\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_wcs = wcs.WCS(naxis=2)\n",
    "my_wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']\n",
    "my_wcs.wcs.crpix = [512, 512]\n",
    "my_wcs.wcs.crval = [70., 20.]\n",
    "my_wcs.wcs.cdelt = [-1/3600, 1/3600]\n",
    "my_wcs.array_shape = [1024, 1024] # NAXIS2, NAXIS1\n",
    "my_wcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting a WCS into an astropy.io.fits.Header instance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The WCS object can be changed and the new WCS can be written out as a new header.\n",
    "\n",
    "By default only the primary WCS keywords are written out to the header. Pass a keyword `relax=True` to write out the SIP distortion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The original WCS\n",
    "w.printwcs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.wcs.crpix = [200, 200]\n",
    "\n",
    "# Calling *to_header()* without arguments writes\n",
    "# out the standard WCS keywords.\n",
    "w.to_header()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Passing *relax=True* writes out the SIP coefficients\n",
    "# and all other distortions.\n",
    "w.to_header(relax=True)"
   ]
  },
  {
   "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(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ref_catalog = Table.read(...)"
   ]
  },
  {
   "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(...)"
   ]
  },
  {
   "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(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gaia_coords = SkyCoord(...)"
   ]
  },
  {
   "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",
    "# Uncomment this block to plot the Gaia coordinates\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",
    "# Repeat for `initial_coords` but use `edgecolors='blue'"
   ]
  },
  {
   "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(...)"
   ]
  },
  {
   "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. Use `sip_degree=3` to fit 3rd-order polynomial distortion.\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(...)"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "### APPENDIX: Optional material on TPV distortion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The TPV distortion convention understood by SExtractor and related programs is also understood by `astropy.wcs`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A header can be read from a text file and used to create the WCS. Here we use a header from an\n",
    "image from the Zwicky Transient Facility, a time-domain imaging survey."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tpv_header = fits.Header.fromtextfile(os.path.join('data', 'ztf_sci_header.txt'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_tpv = wcs.WCS(tpv_header)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The keywords starting with `PV` specify the distortion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_tpv.to_header()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the TPV case, `all_pix2world` gives the same result as `wcs_pix2world`.\n",
    "\n",
    "This is because the underlying WCSLIB library is handling the distortion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(w_tpv.all_pix2world([[500,1000]], 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(w_tpv.wcs_pix2world([[500,1000]], 1))"
   ]
  }
 ],
 "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
}