03 - World Coordinate Systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# World Coordinate Systems in Astropy\n",
    "\n",
    "__Author: Kaustubh Vaghmare__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you open an image is ds9 or any other image viewer, you will notice that as the mouse pointer moves across various regions in the image, a part of the window will reflect the current coordinates in some coordinate system. This is possible because the image header contains a whole bunch of keywords which can be used to compute the coordinates as a function of the pixel coordinates.\n",
    "\n",
    "What one needs is the coordinate for one pixel and then variation in coordinates for every pixel of movement. But if the image is very large, the curvature of the sky can become important. Then additional information to account for second order effects will need to be recorded\n",
    "\n",
    "Let's see how a good fraction of all this is handled using Astropy.\n",
    "\n",
    "Simple steps:\n",
    "\n",
    "    Get hold of the header.\n",
    "    Pass it to a WCS object constructor.\n",
    "    And use this object for coordinate conversions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WCS Keywords\n",
      "\n",
      "Number of WCS axes: 2\n",
      "CTYPE : 'RA---TAN'  'DEC--TAN'  \n",
      "CRVAL : 205.4864335  47.26686298333333  \n",
      "CRPIX : -3929.556588388735  377.0604456526553  \n",
      "PC1_1 PC1_2  : 0.02528041441746252  2.266833305288208e-05  \n",
      "PC2_1 PC2_2  : -1.983416356621748e-05  0.02528846844193638  \n",
      "CDELT : -0.01868178756738699  0.01868178756738699  \n",
      "NAXIS : 777  777\n"
     ]
    }
   ],
   "source": [
    "from astropy.io import fits\n",
    "import astropy.wcs as wcs\n",
    "\n",
    "hdulist = fits.open('data/example2.fits')\n",
    "w = wcs.WCS(hdulist[0].header)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array(202.69037279), array(47.10039636)]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Convert pixels to coordinates.\n",
    "w.wcs_pix2world(100, 100, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([202.69037279, 202.68677899]), array([47.10039636, 47.10266995])]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.wcs_pix2world([100, 105], [100, 105], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array(100.56778084), array(99.1810579)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.wcs_world2pix(202.69, 47.1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[202.76146626,  47.05535561],\n",
       "       [202.74206475,  47.42164629],\n",
       "       [202.2009436 ,  47.40718098],\n",
       "       [202.22406523,  47.04102483]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.calc_footprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some convenience functions\n",
    "\n",
    "There are some utilities provided for your convenience which make coordinate to pixel and vice versa calculations much simple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
       "    (202.68965405, 47.10085109)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output = wcs.utils.pixel_to_skycoord(100, 100, w)\n",
    "output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember: This is capable of ignoring higher dimensions in WCS in radio data, for example.\n",
    "\n",
    "The output here is a special object known as a sky coordinate object. This can be created by saying,\n",
    "\n",
    "    from astropy.coordinates import SkyCoord\n",
    "    c = SkyCoord(10, 20, unit='deg')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<SkyCoord (Galactic): (l, b) in deg\n",
      "    (104.37183205, 68.58518637)>\n"
     ]
    }
   ],
   "source": [
    "import astropy\n",
    "\n",
    "new = output.transform_to(astropy.coordinates.Galactic)\n",
    "print(new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0d28m50.1524s\n"
     ]
    }
   ],
   "source": [
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "distance = output.separation(SkyCoord(202, 47, unit='deg'))\n",
    "print(distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "202d41m22.7546s 47d06m03.0639s\n"
     ]
    }
   ],
   "source": [
    "print(output.ra, output.dec)"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}