04-Coordinates/astropy_coordinates_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An overview of `astropy.coordinates`:\n",
    "\n",
    "### Matching catalogs, planning observations, and handling stellar velocities\n",
    "\n",
    "### *Note: This notebook is a copy of the tutorial notebook with some redundant cells omitted and with exercise solutions filled in*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preliminary imports\n",
    "\n",
    "We start by importing some general packages we will need below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: If the next cell fails, make sure you have installed `astroquery` by running:\n",
    "```\n",
    "pip install astroquery\n",
    "```\n",
    "in your terminal. If so, you will need to go to \"Kernel->Restart\" in this notebook after installing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from astropy.io import fits\n",
    "from astropy.table import QTable\n",
    "from astropy.utils.data import download_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Representing on-sky locations with `astropy.coordinates`\n",
    "\n",
    "In Astropy, the most common way of representing and working with sky coordinates is to use the `SkyCoord` object. A `SkyCoord` can be created directly from angles or arrays of angles with associated units, as demonstrated below. \n",
    "\n",
    "To get started, let's assume that we want to create a `SkyCoord` object for the center of the open cluster NGC 188 so that later we can query and retrieve stars that might be members of the cluster. Let's also assume, for now, that we already know the sky coordinates of the cluster to be (12.11, 85.26) degrees in the ICRS coordinate frame. The ICRS — sometimes referred to simply as \"equitorial\" or \"J2000\" coordinates — is currently the most common astronomical coordinate frame for stellar or extragalactic astronomy, and is the default coordinate frame for `SkyCoord`. We could therefore create a `SkyCoord` object for the center of NGC 188 as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center = SkyCoord(12.11*u.deg, 85.26*u.deg)\n",
    "ngc188_center"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is generally recommended to explicitly specify the frame your coordinates are in; In this case, because ICRS is the default, this would be an equivalent way of creating our `SkyCoord` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center = SkyCoord(12.11*u.deg, 85.26*u.deg, frame='icrs')\n",
    "ngc188_center"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we will see below, there are many other supported coordinate frames.\n",
    "\n",
    "In the above initializations, we passed in `astropy.units.Quantity` objects with angular units to specify the angular components of our sky coordinates. However, `SkyCoord` will also accept string-formatted coordinates either as separate strings for RA/Dec or a single string:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SkyCoord('00h48m26.4s', '85d15m36s', frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some string representations do not explicitly define units, so you would have to specify the units explicitly if they are not included in the string:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SkyCoord('00:48:26.4 85:15:36', unit=(u.hour, u.deg), \n",
    "         frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information and examples on initializing `SkyCoord` objects, [see this documentation](http://docs.astropy.org/en/latest/coordinates/skycoord.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `SkyCoord` object we created has various ways of accessing the information contained within it. All `SkyCoord` objects have attributes that allow you to retrieve the coordinate component data, but the component names will change depending on the coordinate frame of the `SkyCoord` you have. In our case, we have a `SkyCoord` with ICRS coordinates: The component names are therefore lower-case abbreviations of Right Ascension, `.ra`, and Declination, `.dec`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.ra, ngc188_center.dec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These ``ra`` and ``dec`` attributes return specialized [``Quantity``](http://docs.astropy.org/en/stable/units/index.html) objects (actually, a subclass called [``Angle``](http://docs.astropy.org/en/stable/api/astropy.coordinates.Angle.html), which in turn is subclassed by [``Latitude``](http://docs.astropy.org/en/stable/api/astropy.coordinates.Latitude.html) and [``Longitude``](http://docs.astropy.org/en/stable/api/astropy.coordinates.Longitude.html)).  These objects store angles and provide pretty representations of those angles, as well as some useful attributes to quickly convert to common angle units. For example, in a Jupyter notebook, these objects know how to represent themselves using Latex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.ra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.dec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also retrieve the component values in different units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.ra.hour, ngc188_center.ra.radian, ngc188_center.ra.degree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And format the values into strings with specified units ([docs](http://docs.astropy.org/en/latest/coordinates/formatting.html)), e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.ra.to_string(unit=u.hourangle, sep=':')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a `SkyCoord` object for the center of NGC 188, we can select sources from the *Gaia* Data Release 2 catalog around this position to try to find stars that might be members of the cluster. To do this, we will use the `astroquery.gaia` module to query the *Gaia* data archive ([docs](https://astroquery.readthedocs.io/en/latest/gaia/gaia.html)).\n",
    "\n",
    "This requires an internet connection, but if it fails, the catalog file is included in the repository so you can just load it locally (i.e. skip the next cell if you do not have an internet connection):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.gaia import Gaia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = QTable.read('gaia_results.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This table now contains about 5000 stars from *Gaia* DR2 around the coordinate position of the center of NGC 188. Let's now construct a `SkyCoord` object with the results table. In the *Gaia* data archive, the ICRS coordinates of a source are given as column names `\"ra\"` and `\"dec\"`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table['ra']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table['dec']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, because the *Gaia* archive provides data tables with associated units, and we read this table with the `QTable` object ([docs](http://docs.astropy.org/en/latest/table/mixin_columns.html#quantity-and-qtable)), the above columns are represented as `Quantity` objects with units of degrees. Note also that these columns contain many (>5000!) coordinate values: We can pass these directly in to `SkyCoord` to get a single `SkyCoord` object to represent all of these coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_coords = SkyCoord(table['ra'], table['dec'])\n",
    "gaia_coords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a `SkyCoord` of the open cluster The Pleiades (either by looking up the coordinates and passing them in, or by using the convenience method we learned about above):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pleiades_center = SkyCoord.from_name('Pleiades')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using only a single method/function call on the `SkyCoord` object representing the center of NGC 188, print a string with the RA/Dec in the form 'HH:MM:SS.S DD:MM:SS.S'. Check your answer against [SIMBAD](http://simbad.u-strasbg.fr/simbad/), which will show you sexigesimal coordinates for the object.\n",
    "\n",
    "<sub>(Hint: `SkyCoord.to_string()` might be useful)</sub>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.to_string(style=\"hmsdms\", sep=\":\", precision=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a single method/function call on the `SkyCoord` object containing the results of our *Gaia* query, compute the angular separation between each resulting star and the coordinates of the cluster center.\n",
    "\n",
    "<sub>(Hint: `SkyCoord.separation()` might be useful)</sub>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_coords.separation(ngc188_center)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## More than just sky positions: Including distance information in `SkyCoord`\n",
    "\n",
    "So far, we have only used `SkyCoord` to represent sky positions (i.e. `ra` and `dec` only). It is sometimes useful to include distance information with the sky coordinates, thereby fully specifying the 3D position of a given source. To pass in distance information, `SkyCoord` accepts a keyword argument `distance`. So, if we also knew the distance to NGC 188, we could also pass in a distance (as a `Quantity` object) using this argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center_3d = SkyCoord(12.11*u.deg, 85.26*u.deg,\n",
    "                            distance=1.96*u.kpc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our results from querying the *Gaia* catalog, we are also given parallax measurements to each star in the catalog. The default way of passing in a distance, as above, is to pass in a `Quantity` with a unit of length. However, `astropy.coordinates` also provides a specialized object, `Distance`, for handling common transformations of different distance representations ([docs](http://docs.astropy.org/en/latest/coordinates/index.html#distance)). Among other things, this class supports passing in a parallax value <sub>(but note that the transformation to distance here just inverts the parallax, which is not always a good idea!)<sub>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Distance(parallax=1*u.mas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The catalog of stars we queried from *Gaia* also contains parallax information, so we can also create a `SkyCoord` object to represent the 3D positions of all of the *Gaia* stars. Here, we sub-select the table to return only stars with positive parallax values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_3d = table[table['parallax'] > 0]\n",
    "\n",
    "gaia_dist = Distance(parallax=table_3d['parallax'])\n",
    "gaia_coords_3d = SkyCoord(table_3d['ra'], table_3d['dec'],\n",
    "                          distance=gaia_dist)\n",
    "gaia_coords_3d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have 3D position information for both the cluster center, and for the stars we queried from *Gaia*, we can compute the 3D separation (distance) between all of the *Gaia* sources and the cluster center:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sep3d = gaia_coords_3d.separation_3d(ngc188_center_3d)\n",
    "sep3d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "Using the 3D separation values, define a boolean mask to select candidate members of the cluster. Select all stars within 50 pc of the cluster center. How many candidate members of NGC 188 do we have, based on their 3D positions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_3d_mask = sep3d < 50*u.pc\n",
    "ngc188_3d_mask.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Cross-matching and comparing catalogs\n",
    "\n",
    "Let's recap what we have done so far: We defined a `SkyCoord` object to represent the center of an open cluster (NGC 188), we queried the *Gaia* DR2 catalog to select stars that are close (on the sky) to the center of the cluster, and we used the parallax values from *Gaia* to select stars that are near NGC 188 in 3D position.\n",
    "\n",
    "Imagine that we now want to cross-match our catalog of candidate members of NGC 188 to another catalog: Here we will cross-match the *Gaia* stars with the 2MASS photometric catalog to retrieve infrared magnitudes for these stars, and then plot a color-magnitude diagram. To do this, we first need to query the 2MASS catalog to retrieve all sources in a region around the center of NGC 188, as we did for *Gaia*. We will once again use `astroquery` to execute this query: This will again require an internet connection, but we have again included the results of this query in a file along with this notebook in case you are not connected to the internet.\n",
    "\n",
    "To query 2MASS, we will use the `astroquery.vizier` module ([docs](https://astroquery.readthedocs.io/en/latest/vizier/vizier.html)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.vizier import Vizier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmass_table = QTable.read('2MASS_results.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with the *Gaia* results table, we can now create a single `SkyCoord` object to represent all of the sources returned from our query to the 2MASS catalog. Let's look at the column names in this table by displaying the first few rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmass_table[:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From looking at the column names, it looks like the two relevant sky coordinate columns are `RAJ2000` for `ra` and `DEJ2000` for `dec`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmass_coords = SkyCoord(tmass_table['RAJ2000'], \n",
    "                        tmass_table['DEJ2000'])\n",
    "len(tmass_coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to cross-match our *Gaia*-selected candidate members of NGC 188 with this table of photometry from 2MASS. Let's first use the mask we defined above to select members of NGC 188 based on 3D position:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_ngc188_members = gaia_coords_3d[ngc188_3d_mask]\n",
    "gaia_ngc188_members_table = table_3d[ngc188_3d_mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the ``SkyCoord.match_to_catalog_sky`` method to match the two catalogs ([docs](http://docs.astropy.org/en/latest/coordinates/matchsep.html#astropy-coordinates-matching)). Note that order matters: Here we are matching 2MASS to *Gaia* because we ultimately want to make a color-magnitude diagram using the 2MASS photometry. `SkyCoord.match_to_catalog_sky` returns three objects: first, the indices into `gaia_ngc188_members` that get the closest matches in tmass_coords, second, the angular separation between each `tmass_coords` coordinate and the closest source in `gaia_ngc188_members`, and finally, the 3D distance between each `tmass_coords` coordinate and the closest source in `gaia_ngc188_members`. Here, the 3D distances will not be useful because the 2MASS coordinates do not have associated distance information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_tmass, d2d_tmass, d3d_tmass = tmass_coords.match_to_catalog_sky(gaia_ngc188_members)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at the distribution of separations for all of the cross-matched sources:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(d2d_tmass.arcsec, histtype='step', \n",
    "         bins=np.logspace(-2, 2., 64))\n",
    "plt.xlabel('separation [arcsec]')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this, it looks like the sources that actually match to our *Gaia* NGC 188 members are all within about 2 arcseconds, so that is promising! How many of the 2MASS sources matched to one of our NGC 188 candidate members?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(d2d_tmass < 2*u.arcsec).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is a large fraction of the 321 members we identified. Missing sources may be too faint for 2MASS to detect, so we do not expect all of the *Gaia* sources to be in 2MASS. With our cross-match done, we can now make 2MASS+Gaia color-magnitude diagrams of our candidate NGC 188 members:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Jmag = tmass_table['Jmag'][d2d_tmass < 2*u.arcsec]\n",
    "Gmag = gaia_ngc188_members_table['phot_g_mean_mag'][idx_tmass[d2d_tmass < 2*u.arcsec]]\n",
    "Bmag = gaia_ngc188_members_table['phot_bp_mean_mag'][idx_tmass[d2d_tmass < 2*u.arcsec]]\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.scatter(Gmag - Jmag, Gmag, \n",
    "           marker='o', color='k', \n",
    "           linewidth=0, alpha=0.5)\n",
    "ax.set_xlabel('$G - J$')\n",
    "ax.set_ylabel('$G$')\n",
    "ax.set_xlim(0, 3)\n",
    "ax.set_ylim(19, 10) # backwards because magnitudes!\n",
    "\n",
    "ax = axes[1]\n",
    "ax.scatter(Bmag - Gmag, Jmag, \n",
    "           marker='o', color='k', \n",
    "           linewidth=0, alpha=0.5)\n",
    "ax.set_xlabel('$G_{BP} - G$')\n",
    "ax.set_ylabel('$J$')\n",
    "ax.set_xlim(0.2, 1)\n",
    "ax.set_ylim(17, 8) # backwards because magnitudes!\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those both look like color-magnitude diagrams of a main sequence + red giant branch of an intermediate-age stellar cluster, so it looks like our selection and cross-matching has worked!\n",
    "\n",
    "For more on what matching options are available, check out the [separation and matching section of the astropy documentation](https://astropy.readthedocs.io/en/stable/coordinates/matchsep.html).  Or for more on what you can do with `SkyCoord`, see [its API documentation](http://astropy.readthedocs.org/en/stable/api/astropy.coordinates.SkyCoord.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: skip this exercise if you do not have an internet connection*\n",
    "\n",
    "Using the coordinates of the center of the Pleiades you defined in the previous set of exercises, download a catalog of 2MASS sources within 1 degree of the center of the Pleiades. Make a color-magnitude diagram of J-H vs. J for all sources. Can you see the main sequence corresponding to the Pleiades?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = Vizier(catalog=\"II/246\")  \n",
    "v.ROW_LIMIT = -1\n",
    "\n",
    "result = v.query_region(pleiades_center, radius=1*u.deg)\n",
    "tmass_table_pleiades = result[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Jmag_p = tmass_table_pleiades['Jmag']\n",
    "Hmag_p = tmass_table_pleiades['Hmag']\n",
    "\n",
    "plt.scatter(Jmag_p - Hmag_p, Jmag_p,\n",
    "            marker='o', color='k', \n",
    "            linewidth=0, alpha=0.5)\n",
    "plt.xlim(-0.4, 0.75)\n",
    "plt.ylim(12, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: skip this exercise if you do not have an internet connection*\n",
    "\n",
    "Again using the coordinates of the center of the Pleiades, now download all *Gaia* sources within 1 degree of the center of the Pleiades. Cross-match the 2MASS and *Gaia* sources, and make a G-J vs. G color-magnitude diagram of all sources within 200 parsecs of the Sun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = Gaia.cone_search_async(pleiades_center, radius=1*u.deg)\n",
    "gaia_table_pleiades = job.get_results()\n",
    "\n",
    "# This might produce a bunch of warnings: you can ignore these!\n",
    "gaia_table_pleiades = QTable(gaia_table_pleiades[gaia_table_pleiades['parallax'] > 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmass_ple_coords = SkyCoord(tmass_table_pleiades['RAJ2000'],\n",
    "                            tmass_table_pleiades['DEJ2000'])\n",
    "\n",
    "gaia_ple_coords = SkyCoord(gaia_table_pleiades['ra'],\n",
    "                           gaia_table_pleiades['dec'],\n",
    "                           Distance(parallax=gaia_table_pleiades['parallax']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_ple, d2d_ple, _ = gaia_ple_coords.match_to_catalog_sky(tmass_ple_coords)\n",
    "\n",
    "sep_mask = (d2d_ple < 2*u.arcsec) & (gaia_ple_coords.distance < 200*u.pc)\n",
    "\n",
    "Gmag_ple = gaia_table_pleiades[sep_mask]['phot_g_mean_mag']\n",
    "Jmag_ple = tmass_table_pleiades[idx_ple[sep_mask]]['Jmag']\n",
    "\n",
    "plt.scatter(Gmag_ple - Jmag_ple, Gmag_ple, \n",
    "            marker='o', color='k', \n",
    "            linewidth=0, alpha=0.5)\n",
    "plt.xlabel('$G - J$')\n",
    "plt.ylabel('$G$')\n",
    "plt.xlim(-1, 3.5)\n",
    "plt.ylim(19, 0) # backwards because magnitudes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Transforming between coordinate systems and planning observations\n",
    "\n",
    "Imagine that we found something exciting in our exploration of the data for NGC 188, and we would like to follow-up with our own observations of a given sources. `astropy.coordinates` provides tools to assist in planning observations, especially from ground-based observatories, as we will see below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducting frame transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand the code in this section, it may help to read over the [overview of the astropy coordinates scheme](http://astropy.readthedocs.org/en/stable/coordinates/index.html#overview-of-astropy-coordinates-concepts).  The key concept to understand is that all coordinates in astropy are in particular reference \"frames\", and we can transform between a specific `SkyCoord` object in one frame to another. For example, we can transform our previously-defined center of NGC 188 from ICRS to another coordinate frame, like Galactic coordinates. In `astropy.coordinates`, this is typically done using the `.transform_to()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.transform_to('galactic')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method either takes a string name of a coordinate frame, or an instance of a coordinate frame object itself. For example, the class representing the Galactic frame is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import Galactic\n",
    "ngc188_center.transform_to(Galactic())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For coordinate frames that do not require additional parameters, the transformations can also be accessed by attributes of the `SkyCoord` object. Another way to transform to the Galactic frame, then, is to do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.galactic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that changing frames also changes some of the attributes of the object, most often the names of the component attributes. Recall that for the ICRS frame the components were `ra` and `dec`; For the Galactic frame, the corresponding longitude and latitude are named `l` and `b`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_gal = ngc188_center.galactic\n",
    "ngc188_gal.l, ngc188_gal.b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the names `.ra` and `.dec` will therefore *not* work on this transformed `SkyCoord`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# THIS SHOULD ERROR:\n",
    "ngc188_gal.ra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using frame transformations to get to AltAz\n",
    "\n",
    "To determine whether a target is observable from a given observatory on Earth, we need to convert our coordinate to a frame local to an on-earth observer. The most common choice for such a frame is horizontal coordinates, or \"AltAz\" coordinates. The `AltAz` coordinate frame requires a few extra pieces of information to define: a location on Earth, and the time of observation. These are handled with the `EarthLocation` and `Time` objects, respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import EarthLocation\n",
    "from astropy.time import Time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assume we will be observing at Kitt Peak in a few weeks. We can easily retrieve the Earth location of Kitt Peak using the `EarthLocation.of_site` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observing_location = EarthLocation.of_site('Kitt Peak')\n",
    "# List of all locations can be viewed using\n",
    "# EarthLocation.get_site_names()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's pretend our observing run starts on Jan 20 this year:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observing_time = Time('2020-01-20 1:00')  # 1am UTC=6pm AZ mountain time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use these to create an `AltAz` frame object. Note that this frame has even more information about the atmosphere, which can be used to correct for atmospheric refraction. Here we leave that alone, because the default is to ignore this effect (by setting the air pressure to 0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import AltAz\n",
    "\n",
    "altaz = AltAz(location=observing_location, obstime=observing_time)\n",
    "altaz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can just transform our ICRS `SkyCoord` to `AltAz` to get the location in the sky over Kitt Peak at the requested time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.transform_to(altaz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright, the cluster NGC 188 is up at 6PM local time at Kitt Peak, but that is pretty early to be observing. We could just try various times one at a time to see if the airmass is at a darker time, but we can do better: We can make an airmass plot for the whole night:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this gives a Time object with an *array* of times\n",
    "delta_hours = np.linspace(0, 6, 100)*u.hour\n",
    "full_night_times = observing_time + delta_hours\n",
    "full_night_altaz = AltAz(location=observing_location, obstime=full_night_times)\n",
    "full_night_altaz_coords = ngc188_center.transform_to(full_night_altaz)\n",
    "\n",
    "plt.plot(delta_hours, full_night_altaz_coords.secz)\n",
    "plt.xlabel('Hours from 6PM AZ time')\n",
    "plt.ylabel('Airmass [Sec(z)]')\n",
    "plt.ylim(0.9, 3)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like our source is observable, but at relatively high airmass, for most of the night, but the minimum airmass is at the start of the night. When can we start observing with proper dark skies?  Astropy also provides a ``get_sun`` function that can be used to check the sky position of the sun to check, e.g., when we are in 18-degree twilight:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import get_sun\n",
    "\n",
    "full_night_sun_coos = get_sun(full_night_times).transform_to(full_night_altaz)\n",
    "plt.plot(delta_hours, full_night_sun_coos.alt.deg)\n",
    "plt.axhline(-18, color='k')\n",
    "plt.xlabel('Hours from 6PM AZ time')\n",
    "plt.ylabel('Sun altitude')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like it's just below -18 degrees a little after 7PM local time, so we can start observing our target then!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try converting the NGC 188 center coordinates to an equatorial frame at some other equinox a while in the past (e.g., J1900). Do you see the precession of the equinoxes?\n",
    "\n",
    "<sub>Hint: To see a diagram of the supported frames, see [here](http://docs.astropy.org/en/stable/coordinates/#module-astropy.coordinates). `FK5` has a parameter called 'equinox' - that looks promising!</sub>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import FK5\n",
    "\n",
    "ngc188_center.transform_to(FK5(equinox='J2000'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_center.transform_to(FK5(equinox='J1900'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What time is Sunset on our observing night? Try solving for the time at which the Sun's altitude is closest to 0 (i.e. do not just eyeball the plot, like we did above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import Angle\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "def function_to_minimize(local_time):\n",
    "    time = observing_time + local_time*u.hour\n",
    "    aa_frame = AltAz(location=observing_location, obstime=time)\n",
    "    coo = get_sun(time).transform_to(aa_frame)\n",
    "    return coo.alt.deg**2\n",
    "\n",
    "optimization_result = minimize(function_to_minimize, [1]) # Initial guess: 1 hour after 6pm\n",
    "\n",
    "# this gives the local arizona time in PM, because the function we optimized uses the offset from 6pm\n",
    "sunset_time_dec = 6. + optimization_result.x[0]  \n",
    "Angle(sunset_time_dec*u.hour)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## More than just sky positions: Including velocity information in `SkyCoord`\n",
    "\n",
    "As we have seen above, the `SkyCoord` object can store scalars or arrays of positional coordinate information, can transform between different coordinate frames, and can cross-match sets of coordinates. But wait, there's more! `astropy.coordinates` also supports representing and transforming *velocity* information along with positions ([docs](http://docs.astropy.org/en/latest/coordinates/velocities.html)).\n",
    "\n",
    "### Passing velocity data in to `SkyCoord`\n",
    "\n",
    "Velocity data is passed in to `SkyCoord` in a very similar way to the positional coordinates demonstrated above. However, velocity data must be passed in using keyword arguments, meaning that you must write out the velocity component names, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SkyCoord(ra=10*u.deg, \n",
    "         dec=20*u.deg,\n",
    "         pm_ra_cosdec=1*u.mas/u.yr,\n",
    "         pm_dec=2*u.mas/u.yr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, we create a `SkyCoord` instance with sky position and proper motion data. The proper motion component names for ICRS coordinates are, by default, `pm_ra_cosdec` and `pm_dec` (\"pm\" for \"proper motion\"). You can also specify radial velocity data via the `radial_velocity` argument, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_coord = SkyCoord(ra=10*u.deg, \n",
    "                          dec=20*u.deg,\n",
    "                          pm_ra_cosdec=1*u.mas/u.yr,\n",
    "                          pm_dec=2*u.mas/u.yr,\n",
    "                          radial_velocity=100*u.km/u.s)\n",
    "velocity_coord"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `SkyCoord` object with velocity data can be transformed to other frames just like the objects we are used to from above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_coord.transform_to('galactic')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important caveat, however, is that some reference frames require knowing distance in order to transform all velocity components - something we did not specify in the `velocity_coord` definition above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# THIS SHOULD ERROR:\n",
    "velocity_coord.transform_to(altaz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evolving coordinate positions between epochs\n",
    "\n",
    "For nearby or fast-moving stars, a star's position could change appreciably between two well-spaced observations of the source. For such cases, it might be necessary to evolve the position of the star using the proper motion or velocity of the star in order to accurately predict its position at a new time. Below, we will demonstrate this using data from the Digitized Sky Survey (DSS; digital scans of photographic plates observed in the 1950s) and data from the *Gaia* mission (at epoch J2015.5).\n",
    "\n",
    "The star HD 219829 has a very large proper motion, close to 0.5 arcsec/year. Between the DSS and *Gaia*, we therefore expect that the position of the star has changed by about 0.5 arcmin; Let's see if this is the case! To start, we will query the *Gaia* catalog to retrieve data for this star (but skip the cell below if you do not have an internet connection):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hd219829_coord = SkyCoord(ra=349.72896716*u.deg, \n",
    "                          dec=5.40511585*u.deg,\n",
    "                          distance=34.47896*u.pc,\n",
    "                          pm_ra_cosdec=483.41659*u.mas/u.yr,\n",
    "                          pm_dec=-114.86339*u.mas/u.yr,\n",
    "                          obstime=Time('J2015.5'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note above that we specify a new keyword argument, `obstime`, to `SkyCoord`: This is the observation time of the coordinates, and will be important later when we evolve the position of this star using its proper motion.\n",
    "\n",
    "We now have the position and proper motion for HD 219829. Let's now query the DSS to retrieve a FITS image of the field around this star, using the STSCI DSS image cutout service. Again, skip the cell below if you do not have an internet connection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dss_cutout_filename = 'dss_hd219829.fits'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now load the FITS image of the cutout and use `astropy.visualization` to display the image using its World Coordinate System (WCS) info ([docs](http://docs.astropy.org/en/latest/visualization/wcsaxes/index.html)). By passing in the WCS information (included in the FITS cutout header), we can over-plot a marker for the *Gaia*-measured sky position of HD 219829:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.wcs import WCS\n",
    "\n",
    "hdu = fits.open(dss_cutout_filename)[0]\n",
    "wcs = WCS(hdu.header)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8), \n",
    "                       subplot_kw=dict(projection=wcs))\n",
    "ax.imshow(hdu.data, origin='lower', cmap='Greys_r')\n",
    "ax.set_xlabel('RA')\n",
    "ax.set_ylabel('Dec')\n",
    "ax.set_autoscale_on(False)\n",
    "\n",
    "ax.scatter(hd219829_coord.ra.degree,\n",
    "           hd219829_coord.dec.degree,\n",
    "           s=500,\n",
    "           transform=ax.get_transform('world'),\n",
    "           facecolor='none', linewidth=2, color='tab:red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bright star (as observed by DSS) is our target, and the red circle is where *Gaia* observed this star: It has moved! But, we can account for this, and predict the position of the star at around the time the DSS plate was observed (let's assume it was exactly 1950 - this is not strictly correct, but should be good enough).\n",
    "\n",
    "To account for the proper motion of the source and evolve the position to a new time, we can use the `SkyCoord.apply_space_motion()` method ([docs](http://docs.astropy.org/en/latest/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord.apply_space_motion)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hd219829_coord_1950 = hd219829_coord.apply_space_motion(new_obstime=Time('J1950'))  # ignore the ErfaWarning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8), \n",
    "                       subplot_kw=dict(projection=wcs))\n",
    "ax.imshow(hdu.data, origin='lower', cmap='Greys_r')\n",
    "ax.set_xlabel('RA')\n",
    "ax.set_ylabel('Dec')\n",
    "ax.set_autoscale_on(False)\n",
    "\n",
    "ax.scatter(hd219829_coord.ra.degree,\n",
    "           hd219829_coord.dec.degree,\n",
    "           s=500,\n",
    "           transform=ax.get_transform('world'),\n",
    "           facecolor='none', linewidth=2, color='tab:red')\n",
    "\n",
    "# Plot the predicted (past) position:\n",
    "ax.scatter(hd219829_coord_1950.ra.degree,\n",
    "           hd219829_coord_1950.dec.degree,\n",
    "           s=500,\n",
    "           transform=ax.get_transform('world'),\n",
    "           facecolor='none', linewidth=2, color='tab:blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The predicted position of the star from 1950 is much closer to where the star is in the DSS image: Success!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *Gaia* data table we downloaded for sources around NGC 188 also contains proper motion information. Do any of the sources in the field around NGC 188 have large enough proper motions that we might detect a difference between the DSS position of a star and the *Gaia* position? How many stars have an apparent change in position > 5 arcsec over 60 years?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc188_pm = np.sqrt(table['pmra']**2 + table['pmdec']**2)\n",
    "((ngc188_pm * 60*u.year) > 5*u.arcsec).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download a DSS image for the region around the fastest-moving star in the NGC 188 field. By accounting for the proper motion, can you predict where the star is in the DSS image using only the *Gaia* position and proper motion?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrap-up\n",
    "\n",
    "This tutorial covers a lot of material, but `astropy.coordinates` has even more functionality that we were unable to cover in this workshop. For documentation on other features of `astropy.coordinates`, check out [the astropy.coordinates section of the Astropy documentation](http://astropy.readthedocs.org/en/stable/coordinates/index.html).\n",
    "\n",
    "You might also be interested in [the astroplan affiliated package](http://astroplan.readthedocs.org/), which uses the `astropy.coordinates` to do more advanced versions of the observation planning tasks briefly covered above."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "astropy-tutorials": {
   "author": "Erik Tollerud <erik.tollerud@gmail.com>",
   "date": "July 2015",
   "description": "Demonstrates use of astropy.coordinates for common tasks. Includes matching catalogs against each other, basic observing planning tasks, and basic usage of coordinates.",
   "link_name": "Using astropy.coordinates to Match Catalogs and Plan Observations",
   "name": "",
   "published": true
  },
  "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.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}