{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "title: \"Photometry\"\n",
    "teaching: 3000\n",
    "exercises: 0\n",
    "questions:\n",
    "\n",
    "- \"How do we use Matplotlib to select a polygon and Pandas to merge data from multiple tables?\"\n",
    "\n",
    "objectives:\n",
    "\n",
    "- \"Use Matplotlib to specify a polygon and determine which points fall inside it.\"\n",
    "\n",
    "- \"Use Pandas to merge data from multiple `DataFrames`, much like a database `JOIN` operation.\"\n",
    "\n",
    "keypoints:\n",
    "\n",
    "- \"Matplotlib provides operations for working with points, polygons, and other geometric entities, so it's not just for making figures.\"\n",
    "\n",
    "- \"If you want to perform something like a database `JOIN` operation with data that is in a Pandas `DataFrame`, you can use the `join` or `merge` function.  In many cases, `merge` is easier to use because the arguments are more like SQL.\"\n",
    "\n",
    "- \"Use Matplotlib options to control the size and aspect ratio of figures to make them easier to interpret.\"\n",
    "\n",
    "- \"Be sure to record every element of the data analysis pipeline that would be needed to replicate the results.\"\n",
    "\n",
    "---\n",
    "\n",
    "{% include links.md %}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Photometry\n",
    "\n",
    "This is the sixth in a series of notebooks related to astronomy data.\n",
    "\n",
    "As a continuing example, we will replicate part of the analysis in a recent paper, \"[Off the beaten path: Gaia reveals GD-1 stars outside of the main stream](https://arxiv.org/abs/1805.00425)\" by Adrian M. Price-Whelan and Ana Bonaca.\n",
    "\n",
    "In the previous lesson we downloaded photometry data from Pan-STARRS, which is available from the same server we've been using to get Gaia data. \n",
    "\n",
    "The next step in the analysis is to select candidate stars based on the photometry data.  \n",
    "The following figure from the paper is a color-magnitude diagram showing the stars we previously selected based on proper motion:\n",
    "\n",
    "<img width=\"300\" src=\"https://github.com/datacarpentry/astronomy-python/raw/gh-pages/fig/gd1-3.png\">\n",
    "\n",
    "In red is a theoretical isochrone, showing where we expect the stars in GD-1 to fall based on the metallicity and age of their original globular cluster. \n",
    "\n",
    "By selecting stars in the shaded area, we can further distinguish the main sequence of GD-1 from mostly younger background stars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    "Here are the steps in this notebook:\n",
    "\n",
    "1. We'll reload the data from the previous notebook and make a color-magnitude diagram.\n",
    "\n",
    "2. We'll use an isochrone computed by MIST to specify a polygonal region in the color-magnitude diagram and select the stars inside it.\n",
    "\n",
    "3. Then we'll merge the photometry data with the list of candidate stars, storing the result in a Pandas `DataFrame`.\n",
    "\n",
    "After completing this lesson, you should be able to\n",
    "\n",
    "* Use Matplotlib to specify a `Polygon` and determine which points fall inside it.\n",
    "\n",
    "* Use Pandas to merge data from multiple `DataFrames`, much like a database `JOIN` operation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "## Installing libraries\n",
    "\n",
    "If you are running this notebook on Colab, you can run the following cell to install the libraries we'll use.\n",
    "\n",
    "If you are running this notebook on your own computer, you might have to install these libraries yourself.  See the instructions in the preface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# If we're running on Colab, install libraries\n",
    "\n",
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "\n",
    "if IN_COLAB:\n",
    "    !pip install wget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reload the data\n",
    "\n",
    "The following cell downloads the photometry data we created in the previous notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from wget import download\n",
    "\n",
    "filename = 'gd1_photo.fits'\n",
    "filepath = 'https://github.com/AllenDowney/AstronomicalData/raw/main/data/'\n",
    "\n",
    "if not os.path.exists(filename):\n",
    "    print(download(filepath+filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can read the data back into an Astropy `Table`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "\n",
    "photo_table = Table.read(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting photometry data\n",
    "\n",
    "Now that we have photometry data from Pan-STARRS, we can replicate the [color-magnitude diagram](https://en.wikipedia.org/wiki/Galaxy_color%E2%80%93magnitude_diagram) from the original paper:\n",
    "\n",
    "<img width=\"300\" src=\"https://github.com/datacarpentry/astronomy-python/raw/gh-pages/fig/gd1-3.png\">\n",
    "\n",
    "The y-axis shows the apparent magnitude of each source with the [g filter](https://en.wikipedia.org/wiki/Photometric_system).\n",
    "\n",
    "The x-axis shows the difference in apparent magnitude between the g and i filters, which indicates color.\n",
    "\n",
    "Stars with lower values of (g-i) are brighter in g-band than in i-band, compared to other stars, which means they are bluer.\n",
    "\n",
    "Stars in the lower-left quadrant of this diagram are less bright than the others, and have lower metallicity, which means they are [likely to be older](http://spiff.rit.edu/classes/ladder/lectures/ordinary_stars/ordinary.html).\n",
    "\n",
    "Since we expect the stars in GD-1 to be older than the background stars, the stars in the lower-left are more likely to be in GD-1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes a table containing photometry data and draws a color-magnitude diagram.\n",
    "The input can be an Astropy `Table` or Pandas `DataFrame`, as long as it has columns named `g_mean_psf_mag` and `i_mean_psf_mag`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_cmd(table):\n",
    "    \"\"\"Plot a color magnitude diagram.\n",
    "    \n",
    "    table: Table or DataFrame with photometry data\n",
    "    \"\"\"\n",
    "    y = table['g_mean_psf_mag']\n",
    "    x = table['g_mean_psf_mag'] - table['i_mean_psf_mag']\n",
    "\n",
    "    plt.plot(x, y, 'ko', markersize=0.3, alpha=0.3)\n",
    "\n",
    "    plt.xlim([0, 1.5])\n",
    "    plt.ylim([14, 22])\n",
    "    plt.gca().invert_yaxis()\n",
    "\n",
    "    plt.ylabel('$Magnitude (g)$')\n",
    "    plt.xlabel('$Color (g-i)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`plot_cmd` uses a new function, `invert_yaxis`, to invert the `y` axis, which is conventional when plotting magnitudes, since lower magnitude indicates higher brightness.\n",
    "\n",
    "`invert_yaxis` is a little different from the other functions we've used.  You can't call it like this:\n",
    "\n",
    "```\n",
    "plt.invert_yaxis()          # doesn't work\n",
    "```\n",
    "\n",
    "You have to call it like this:\n",
    "\n",
    "```\n",
    "plt.gca().invert_yaxis()          # works\n",
    "```\n",
    "\n",
    "`gca` stands for \"get current axis\".  It returns an object that represents the axes of the current figure, and that object provides `invert_yaxis`.\n",
    "\n",
    "**In case anyone asks:** The most likely reason for this inconsistency in the interface is that `invert_yaxis` is a lesser-used function, so it's not made available at the top level of the interface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29124 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_cmd(photo_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our figure does not look exactly like the one in the paper because we are working with a smaller region of the sky, so we don't have as many stars.  But we can see an overdense region in the lower left that contains stars with the photometry we expect for GD-1.\n",
    "\n",
    "In the next section we'll use an isochrone to specify a polygon that contains this overdense regioin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Isochrone\n",
    "\n",
    "Based on our best estimates for the ages of the stars in GD-1 and their metallicity, we can compute a [stellar isochrone](https://en.wikipedia.org/wiki/Stellar_isochrone) that predicts the relationship between their magnitude and color.\n",
    "\n",
    "In fact, we can use [MESA Isochrones & Stellar Tracks](http://waps.cfa.harvard.edu/MIST/) (MIST) to compute it for us.\n",
    "\n",
    "Using the [MIST Version 1.2 web interface](http://waps.cfa.harvard.edu/MIST/interp_isos.html), we computed an isochrone with the following parameters:\n",
    "    \n",
    "* Rotation initial v/v_crit = 0.4\n",
    "\n",
    "* Single age, linear scale = 12e9\n",
    "\n",
    "* Composition [Fe/H] = -1.35\n",
    "\n",
    "* Synthetic Photometry, PanStarrs\n",
    "\n",
    "* Extinction av = 0\n",
    "\n",
    "The following cell downloads the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from wget import download\n",
    "\n",
    "filename = 'MIST_iso_5fd2532653c27.iso.cmd'\n",
    "filepath = 'https://github.com/AllenDowney/AstronomicalData/raw/main/data/'\n",
    "\n",
    "if not os.path.exists(filename):\n",
    "    print(download(filepath+filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To read this file we'll download a Python module [from this repository](https://github.com/jieunchoi/MIST_codes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from wget import download\n",
    "\n",
    "filename = 'read_mist_models.py'\n",
    "filepath = 'https://github.com/jieunchoi/MIST_codes/raw/master/scripts/'\n",
    "\n",
    "if not os.path.exists(filename):\n",
    "    print(download(filepath+filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can read the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading in: MIST_iso_5fd2532653c27.iso.cmd\n"
     ]
    }
   ],
   "source": [
    "import read_mist_models\n",
    "\n",
    "filename = 'MIST_iso_5fd2532653c27.iso.cmd'\n",
    "iso = read_mist_models.ISOCMD(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an `ISOCMD` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "read_mist_models.ISOCMD"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(iso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It contains a list of arrays, one for each isochrone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "list"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(iso.isocmds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We only got one isochrone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(iso.isocmds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we can select it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "iso_array = iso.isocmds[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a NumPy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(iso_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But it's an unusual NumPy array, because it contains names for the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype([('EEP', '<i4'), ('isochrone_age_yr', '<f8'), ('initial_mass', '<f8'), ('star_mass', '<f8'), ('log_Teff', '<f8'), ('log_g', '<f8'), ('log_L', '<f8'), ('[Fe/H]_init', '<f8'), ('[Fe/H]', '<f8'), ('PS_g', '<f8'), ('PS_r', '<f8'), ('PS_i', '<f8'), ('PS_z', '<f8'), ('PS_y', '<f8'), ('PS_w', '<f8'), ('PS_open', '<f8'), ('phase', '<f8')])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iso_array.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which means we can select columns using the bracket operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., ..., 6., 6., 6.])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iso_array['phase']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `phase` to select the part of the isochrone for stars in the main sequence and red giant phases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "354"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phase_mask = (iso_array['phase'] >= 0) & (iso_array['phase'] < 3)\n",
    "phase_mask.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "354"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "main_sequence = iso_array[phase_mask]\n",
    "len(main_sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other two columns we'll use are `PS_g` and `PS_i`, which contain simulated photometry data for stars with the given age and metallicity, based on a model of the Pan-STARRS sensors.\n",
    "\n",
    "We'll use these columns to superimpose the isochrone on the color-magnitude diagram, but first we have to use a [distance modulus](https://en.wikipedia.org/wiki/Distance_modulus) to scale the isochrone based on the estimated distance of GD-1.\n",
    "\n",
    "We can use the `Distance` object from Astropy to compute the distance modulus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.4604730134524"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import astropy.coordinates as coord\n",
    "import astropy.units as u\n",
    "\n",
    "distance = 7.8 * u.kpc\n",
    "distmod = coord.Distance(distance).distmod.value\n",
    "distmod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the scaled magnitude and color of the isochrone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "mag_g = main_sequence['PS_g'] + distmod\n",
    "color_g_i = main_sequence['PS_g'] - main_sequence['PS_i']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot it on the color-magnitude diagram like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35224 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_cmd(photo_table)\n",
    "plt.plot(color_g_i, mag_g);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The theoretical isochrone passes through the overdense region where we expect to find stars in GD-1.\n",
    "\n",
    "Let's save this result so we can reload it later without repeating the steps in this section.\n",
    "\n",
    "So we can save the data in an HDF5 file, we'll put it in a Pandas `DataFrame` first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mag_g</th>\n",
       "      <th>color_g_i</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>28.294743</td>\n",
       "      <td>2.195021</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>28.189718</td>\n",
       "      <td>2.166076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>28.051761</td>\n",
       "      <td>2.129312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>27.916194</td>\n",
       "      <td>2.093721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>27.780024</td>\n",
       "      <td>2.058585</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mag_g  color_g_i\n",
       "0  28.294743   2.195021\n",
       "1  28.189718   2.166076\n",
       "2  28.051761   2.129312\n",
       "3  27.916194   2.093721\n",
       "4  27.780024   2.058585"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "iso_df = pd.DataFrame()\n",
    "iso_df['mag_g'] = mag_g\n",
    "iso_df['color_g_i'] = color_g_i\n",
    "\n",
    "iso_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then save it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'gd1_isochrone.hdf5'\n",
    "\n",
    "iso_df.to_hdf(filename, 'iso_df')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making a polygon\n",
    "\n",
    "The following cell downloads the isochrone we made in the previous section, if necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from wget import download\n",
    "\n",
    "filename = 'gd1_isochrone.hdf5'\n",
    "filepath = 'https://github.com/AllenDowney/AstronomicalData/raw/main/data/'\n",
    "\n",
    "if not os.path.exists(filename):\n",
    "    print(download(filepath+filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can read it back in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mag_g</th>\n",
       "      <th>color_g_i</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>28.294743</td>\n",
       "      <td>2.195021</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>28.189718</td>\n",
       "      <td>2.166076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>28.051761</td>\n",
       "      <td>2.129312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>27.916194</td>\n",
       "      <td>2.093721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>27.780024</td>\n",
       "      <td>2.058585</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mag_g  color_g_i\n",
       "0  28.294743   2.195021\n",
       "1  28.189718   2.166076\n",
       "2  28.051761   2.129312\n",
       "3  27.916194   2.093721\n",
       "4  27.780024   2.058585"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iso_df = pd.read_hdf(filename, 'iso_df')\n",
    "iso_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the isochrone looks like on the color-magnitude diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35224 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_cmd(photo_table)\n",
    "plt.plot(iso_df['color_g_i'], iso_df['mag_g']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the bottom half of the figure, the isochrone passes through the overdense region where the stars are likely to belong to GD-1.\n",
    "\n",
    "In the top half, the isochrone passes through other regions where the stars have higher magnitude and metallicity than we expect for stars in GD-1.\n",
    "\n",
    "So we'll select the part of the isochrone that lies in the overdense region.\n",
    "\n",
    "`g_mask` is a Boolean Series that is `True` where `g` is between 18.0 and 21.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "117"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = iso_df['mag_g']\n",
    "\n",
    "g_mask = (g > 18.0) & (g < 21.5)\n",
    "g_mask.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use it to select the corresponding rows in `iso_df`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mag_g</th>\n",
       "      <th>color_g_i</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>94</th>\n",
       "      <td>21.411746</td>\n",
       "      <td>0.692171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>21.322466</td>\n",
       "      <td>0.670238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>21.233380</td>\n",
       "      <td>0.648449</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>21.144427</td>\n",
       "      <td>0.626924</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>21.054549</td>\n",
       "      <td>0.605461</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mag_g  color_g_i\n",
       "94  21.411746   0.692171\n",
       "95  21.322466   0.670238\n",
       "96  21.233380   0.648449\n",
       "97  21.144427   0.626924\n",
       "98  21.054549   0.605461"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iso_masked = iso_df[g_mask]\n",
    "iso_masked.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to select the stars in the overdense region, we have to define a polygon that includes stars near the isochrone.\n",
    "\n",
    "The original paper uses the following formulas to define the left and right boundaries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = iso_masked['mag_g']\n",
    "left_color = iso_masked['color_g_i'] - 0.4 * (g/28)**5\n",
    "right_color = iso_masked['color_g_i'] + 0.8 * (g/28)**5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The intention is to define a polygon that gets wider as `g` increases, to reflect increasing uncertainty.\n",
    "\n",
    "But we can do about as well with a simpler formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = iso_masked['mag_g']\n",
    "left_color = iso_masked['color_g_i'] - 0.06\n",
    "right_color = iso_masked['color_g_i'] + 0.12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what these boundaries look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35436 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_cmd(photo_table)\n",
    "\n",
    "plt.plot(left_color, g, label='left color')\n",
    "plt.plot(right_color, g, label='right color')\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Which points are in the polygon?\n",
    "\n",
    "Matplotlib provides a `Polygon` object that we can use to check which points fall in the polygon we just constructed.\n",
    "\n",
    "To make a `Polygon`, we need to assemble `g`, `left_color`, and `right_color` into a loop, so the points in `left_color` are connected to the points of `right_color` in reverse order.\n",
    "\n",
    "We'll use the following function, which takes two arrays and joins them front-to-back:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def front_to_back(first, second):\n",
    "    \"\"\"Join two arrays front to back.\"\"\"\n",
    "    return np.append(first, second[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`front_to_back` uses a \"slice index\" to reverse the elements of `second`.\n",
    "\n",
    "As explained in the [NumPy documentation](https://numpy.org/doc/stable/reference/arrays.indexing.html), a slice index has three parts separated by colons:\n",
    "\n",
    "* `start`: The index of the element where the slice starts.\n",
    "\n",
    "* `stop`: The index of the element where the slice ends.\n",
    "\n",
    "* `step`: The step size between elements.\n",
    "\n",
    "In this example, `start` and `stop` are omitted, which means all elements are selected.\n",
    "\n",
    "And `step` is `-1`, which means the elements are in reverse order.\n",
    "\n",
    "We can use `front_to_back` to make a loop that includes the elements of `left_color` and `right_color`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(234,)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "color_loop = front_to_back(left_color, right_color)\n",
    "color_loop.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a corresponding loop with the elements of `g` in forward and reverse order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(234,)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mag_loop = front_to_back(g, g)\n",
    "mag_loop.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the loop looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34012 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_cmd(photo_table)\n",
    "plt.plot(color_loop, mag_loop);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a `Polygon`, it will be convenient to put `color_loop` and `mag_loop` into a `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "loop_df = pd.DataFrame()\n",
    "loop_df['color_loop'] = color_loop\n",
    "loop_df['mag_loop'] = mag_loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can pass `loop_df` to `Polygon`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.patches.Polygon at 0x7fe98cd29400>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from matplotlib.patches import Polygon\n",
    "\n",
    "polygon = Polygon(loop_df)\n",
    "polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a `Polygon` object , which provides `contains_points`, which figures out which points are inside the polygon.\n",
    "\n",
    "To test it, we'll create a list with two points, one inside the polygon and one outside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [(0.4, 20), \n",
    "          (0.4, 16)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can make sure `contains_points` does what we expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True, False])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inside = polygon.contains_points(points)\n",
    "inside"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an array of Boolean values.\n",
    "\n",
    "We are almost ready to select stars whose photometry data falls in this polygon.  But first we need to do some data cleaning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the polygon\n",
    "\n",
    "[Reproducibile research](https://en.wikipedia.org/wiki/Reproducibility#Reproducible_research) is \"the idea that ... the full computational environment used to produce the results in the paper such as the code, data, etc. can be used to reproduce the results and create new work based on the research.\"\n",
    "\n",
    "This Jupyter notebook is an example of reproducible research because it contains all of the code needed to reproduce the results, including the database queries that download the data and and analysis.\n",
    "\n",
    "In this lesson we used an isochrone to derive a polygon, which we used to select stars based on photometry. \n",
    "So it is important to record the polygon as part of the data analysis pipeline.\n",
    "\n",
    "Here's how we can save it in an HDF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'gd1_polygon.hdf5'\n",
    "loop_df.to_hdf(filename, 'loop_df')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reloading the data\n",
    "\n",
    "Now we need to combine the photometry data with the list of candidate stars we identified in a previous notebook.  The following cell downloads it:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from wget import download\n",
    "\n",
    "filename = 'gd1_candidates.hdf5'\n",
    "filepath = 'https://github.com/AllenDowney/AstronomicalData/raw/main/data/'\n",
    "\n",
    "if not os.path.exists(filename):\n",
    "    print(download(filepath+filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "candidate_df = pd.read_hdf(filename, 'candidate_df')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`candidate_df` is the Pandas DataFrame that contains the results from Lesson 4, which selects stars likely to be in GD-1 based on proper motion.  It also includes position and proper motion transformed to the ICRS frame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merging photometry data\n",
    "\n",
    "Before we select stars based on photometry data, we have to solve two problems:\n",
    "\n",
    "1. We only have Pan-STARRS data for some stars in `candidate_df`.\n",
    "\n",
    "2. Even for the stars where we have Pan-STARRS data in `photo_table`, some photometry data is missing.\n",
    "\n",
    "We will solve these problems in two step:\n",
    "\n",
    "1. We'll merge the data from `candidate_df` and `photo_table` into a single Pandas `DataFrame`.\n",
    "\n",
    "2. We'll use Pandas functions to deal with missing data.\n",
    "\n",
    "`candidate_df` is already a `DataFrame`, but `results` is an Astropy `Table`.  Let's convert it to Pandas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "source_id\n",
      "g_mean_psf_mag\n",
      "i_mean_psf_mag\n"
     ]
    }
   ],
   "source": [
    "photo_df = photo_table.to_pandas()\n",
    "\n",
    "for colname in photo_df.columns:\n",
    "    print(colname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to combine `candidate_df` and `photo_df` into a single table, using `source_id` to match up the rows.\n",
    "\n",
    "You might recognize this task; it's the same as the JOIN operation in ADQL/SQL.\n",
    "\n",
    "Pandas provides a function called `merge` that does what we want.  Here's how we use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>source_id</th>\n",
       "      <th>ra</th>\n",
       "      <th>dec</th>\n",
       "      <th>pmra</th>\n",
       "      <th>pmdec</th>\n",
       "      <th>parallax</th>\n",
       "      <th>radial_velocity</th>\n",
       "      <th>phi1</th>\n",
       "      <th>phi2</th>\n",
       "      <th>pm_phi1</th>\n",
       "      <th>pm_phi2</th>\n",
       "      <th>g_mean_psf_mag</th>\n",
       "      <th>i_mean_psf_mag</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>635860218726658176</td>\n",
       "      <td>138.518707</td>\n",
       "      <td>19.092339</td>\n",
       "      <td>-5.941679</td>\n",
       "      <td>-11.346409</td>\n",
       "      <td>0.307456</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-59.247330</td>\n",
       "      <td>-2.016078</td>\n",
       "      <td>-7.527126</td>\n",
       "      <td>1.748779</td>\n",
       "      <td>17.8978</td>\n",
       "      <td>17.517401</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>635674126383965568</td>\n",
       "      <td>138.842874</td>\n",
       "      <td>19.031798</td>\n",
       "      <td>-3.897001</td>\n",
       "      <td>-12.702780</td>\n",
       "      <td>0.779463</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-59.133391</td>\n",
       "      <td>-2.306901</td>\n",
       "      <td>-7.560608</td>\n",
       "      <td>-0.741800</td>\n",
       "      <td>19.2873</td>\n",
       "      <td>17.678101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>635535454774983040</td>\n",
       "      <td>137.837752</td>\n",
       "      <td>18.864007</td>\n",
       "      <td>-4.335041</td>\n",
       "      <td>-14.492309</td>\n",
       "      <td>0.314514</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-59.785300</td>\n",
       "      <td>-1.594569</td>\n",
       "      <td>-9.357536</td>\n",
       "      <td>-1.218492</td>\n",
       "      <td>16.9238</td>\n",
       "      <td>16.478100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>635497276810313600</td>\n",
       "      <td>138.044516</td>\n",
       "      <td>19.009471</td>\n",
       "      <td>-7.172931</td>\n",
       "      <td>-12.291499</td>\n",
       "      <td>0.425404</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-59.557744</td>\n",
       "      <td>-1.682147</td>\n",
       "      <td>-9.000831</td>\n",
       "      <td>2.334407</td>\n",
       "      <td>19.9242</td>\n",
       "      <td>18.334000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>635614168640132864</td>\n",
       "      <td>139.592197</td>\n",
       "      <td>18.807956</td>\n",
       "      <td>-3.309603</td>\n",
       "      <td>-13.708905</td>\n",
       "      <td>0.583382</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-58.938113</td>\n",
       "      <td>-3.024192</td>\n",
       "      <td>-8.062762</td>\n",
       "      <td>-1.869082</td>\n",
       "      <td>16.1516</td>\n",
       "      <td>14.666300</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            source_id          ra        dec      pmra      pmdec  parallax  \\\n",
       "0  635860218726658176  138.518707  19.092339 -5.941679 -11.346409  0.307456   \n",
       "1  635674126383965568  138.842874  19.031798 -3.897001 -12.702780  0.779463   \n",
       "2  635535454774983040  137.837752  18.864007 -4.335041 -14.492309  0.314514   \n",
       "3  635497276810313600  138.044516  19.009471 -7.172931 -12.291499  0.425404   \n",
       "4  635614168640132864  139.592197  18.807956 -3.309603 -13.708905  0.583382   \n",
       "\n",
       "   radial_velocity       phi1      phi2   pm_phi1   pm_phi2  g_mean_psf_mag  \\\n",
       "0              NaN -59.247330 -2.016078 -7.527126  1.748779         17.8978   \n",
       "1              NaN -59.133391 -2.306901 -7.560608 -0.741800         19.2873   \n",
       "2              NaN -59.785300 -1.594569 -9.357536 -1.218492         16.9238   \n",
       "3              NaN -59.557744 -1.682147 -9.000831  2.334407         19.9242   \n",
       "4              NaN -58.938113 -3.024192 -8.062762 -1.869082         16.1516   \n",
       "\n",
       "   i_mean_psf_mag  \n",
       "0       17.517401  \n",
       "1       17.678101  \n",
       "2       16.478100  \n",
       "3       18.334000  \n",
       "4       14.666300  "
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "merged = pd.merge(candidate_df, \n",
    "                  photo_df, \n",
    "                  on='source_id')\n",
    "merged.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first argument is the \"left\" table, the second argument is the \"right\" table, and the keyword argument `on='source_id'` specifies a column to use to match up the rows.\n",
    "\n",
    "The result is a `DataFrame` that contains the same number of rows as `photo_df`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7346, 3724, 3724)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(candidate_df), len(photo_df), len(merged)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And it contains all columns from both tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "source_id\n",
      "ra\n",
      "dec\n",
      "pmra\n",
      "pmdec\n",
      "parallax\n",
      "radial_velocity\n",
      "phi1\n",
      "phi2\n",
      "pm_phi1\n",
      "pm_phi2\n",
      "g_mean_psf_mag\n",
      "i_mean_psf_mag\n"
     ]
    }
   ],
   "source": [
    "for colname in merged.columns:\n",
    "    print(colname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Detail** You might notice that Pandas also provides a function called `join`; it does almost the same thing, but the interface is slightly different.  We think `merge` is a little easier to use, so that's what we chose.  It's also more consistent with JOIN in SQL, so if you learn how to use `pd.merge`, you are also learning how to use SQL JOIN.\n",
    "\n",
    "Also, someone might ask why we have to use Pandas to do this join; why didn't we do it in ADQL.  The answer is that we could have done that, but since we already have the data we need, we should probably do the computation locally rather than make another round trip to the Gaia server."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting based on photometry\n",
    "\n",
    "Now let's see how many of these points are inside the polygon we chose.\n",
    "\n",
    "We'll put color and magnitude data from `merged` into a new `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>color</th>\n",
       "      <th>mag</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.3804</td>\n",
       "      <td>17.8978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.6092</td>\n",
       "      <td>19.2873</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.4457</td>\n",
       "      <td>16.9238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.5902</td>\n",
       "      <td>19.9242</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.4853</td>\n",
       "      <td>16.1516</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    color      mag\n",
       "0  0.3804  17.8978\n",
       "1  1.6092  19.2873\n",
       "2  0.4457  16.9238\n",
       "3  1.5902  19.9242\n",
       "4  1.4853  16.1516"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points = pd.DataFrame()\n",
    "\n",
    "points['color'] = merged['g_mean_psf_mag'] - merged['i_mean_psf_mag']\n",
    "points['mag'] = merged['g_mean_psf_mag']\n",
    "\n",
    "points.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which we can pass to `contains_points`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False, False, False, ..., False, False, False])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inside = polygon.contains_points(points)\n",
    "inside"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a Boolean array.  We can use `sum` to see how many stars fall in the polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "454"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inside.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `inside` as a mask to select stars that fall inside the polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected2 = merged[inside]\n",
    "points2 = points[inside]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a color-magnitude plot one more time, highlighting the selected stars with green markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 42944 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_cmd(photo_table)\n",
    "plt.plot(color_g_i, mag_g)\n",
    "plt.plot(color_loop, mag_loop)\n",
    "\n",
    "plt.plot(points2['color'], points2['mag'], 'g.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the selected stars are, in fact, inside the polygon, which means they have photometry data consistent with GD-1.\n",
    "\n",
    "Finally, we can plot the coordinates of the selected stars:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12916 characters]\n",
      "text/plain": [
       "<Figure size 720x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,2.5))\n",
    "\n",
    "x = selected2['phi1']\n",
    "y = selected2['phi2']\n",
    "\n",
    "plt.plot(x, y, 'ko', markersize=0.7, alpha=0.9)\n",
    "\n",
    "plt.xlabel('ra (degree GD1)')\n",
    "plt.ylabel('dec (degree GD1)')\n",
    "\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example includes two new Matplotlib commands:\n",
    "\n",
    "* `figure` creates the figure.  In previous examples, we didn't have to use this function; the figure was created automatically.  But when we call it explicitly, we can provide arguments like `figsize`, which sets the size of the figure.\n",
    "\n",
    "* `axis` with the parameter `equal` sets up the axes so a unit is the same size along the `x` and `y` axes.\n",
    "\n",
    "In an example like this, where `x` and `y` represent coordinates in space, equal axes ensures that the distance between points is represented accurately.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the data\n",
    "\n",
    "Finally, let's write the merged DataFrame to a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'gd1_merged.hdf5'\n",
    "\n",
    "merged.to_hdf(filename, 'merged')\n",
    "selected2.to_hdf(filename, 'selected2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-rw-r-- 1 downey downey 1.1M Dec 29 11:51 gd1_merged.hdf5\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lh gd1_merged.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are using Windows, `ls` might not work; in that case, try:\n",
    "\n",
    "```\n",
    "!dir gd1_merged.hdf5\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this lesson, we worked with three datasets: \n",
    "\n",
    "* The list of candidate stars from Gaia,\n",
    "\n",
    "* The photometry data from Pan-STARRS, and\n",
    "\n",
    "* An isochrone computed by MIST.\n",
    "\n",
    "We drew a color-magnitude diagram and used it to identify stars we think are likely to be in GD-1.\n",
    "\n",
    "We used the isochrone to define a polygon that includes those stars.\n",
    "\n",
    "Then we used a Pandas `merge` operation to combine Gaia and Pan-STARRS data into a single `DataFrame`.\n",
    "\n",
    "Plotting the results, we have a clear picture of GD-1, similar to Figure 1 in the original paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Best practices\n",
    "\n",
    "* Matplotlib provides operations for working with points, polygons, and other geometric entities, so it's not just for making figures.\n",
    "\n",
    "* If you want to perform something like a database `JOIN` operation with data that is in a Pandas `DataFrame`, you can use the `join` or `merge` function.  In many cases, `merge` is easier to use because the arguments are more like SQL.\n",
    "\n",
    "* Use Matplotlib options to control the size and aspect ratio of figures to make them easier to interpret.  In this example, we scaled the axes so the size of a degree is equal along both axes.\n",
    "\n",
    "* Be sure to record every element of the data analysis pipeline that would be needed to replicate the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
