tutorials/tutorial_02/pysnack_1_02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Astropy Tables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Astropy tables: columns and units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.table import Table, QTable\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 4, 5], dtype=np.int32)\n",
    "b = [2.0, 5.0, 8.5] * u.cm\n",
    "c = ['x', 'y', 'z']\n",
    "d = [10, 20, 30] * u.m / u.s\n",
    "\n",
    "t = QTable([a, b, c, d],\n",
    "           names=('id', 'length', 'label', 'velocity'),\n",
    "           meta={'name': 'first table'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t['velocity']   # Show QTable vs Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t['velocity'][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t['velocity'][0:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t['velocity'].unit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining information and statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t.info('stats')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = t.info('stats', out=None)\n",
    "stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats[3]['mean']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time and coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.time import Time\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "Time.FORMATS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mjd = Time([56200.25, 56400.33, 57500.66, 58000], format='mjd', scale='utc')\n",
    "date = mjd.strftime('%Y-%m-%d %H:%M')   # You can also convert to .to_datetime() object\n",
    "\n",
    "sc = SkyCoord([10, 20, 30, 40], [-45, +40, +55, 33], unit='deg')\n",
    "\n",
    "tab = QTable([mjd, date,  sc, sc.to_string('hmsdms')], names=['MJD', 'Date', 'skycoord', 'coord'])\n",
    "tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['coord']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Managing table columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab.colnames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substitute column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['MJD'] = [56100, 56200, 56300, 56400]\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create new column. Directly or with operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['x'] = [1,2,3,4]\n",
    "tab['y'] = [10,20,30,40]\n",
    "tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['z'] = tab['x']*2 + tab['y']\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Selection and slicing of relevant columns "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['x']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(tab['x'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['x'].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['x'].mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select multiple columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab[['MJD', 'x', 'z']]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_columns = ['x', 'y', 'z']\n",
    "tab[my_columns]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QTable vs Table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⚠️ `Table` with units have limited functionality. For full unit management, it is recommended to use `QTable`. For most applications, `Table` is enough, but if you have to manage and operate with columns with maximum compatibility with unit management, `QTable` is recommended. The caveat is that `QTables` are much more strict and slightly limited when dealing with other python modules, for example to work with matplotlib you constantly need to use `with quantity_support()`. More info in [Plotting Astropy objects in Matplotlib](https://docs.astropy.org/en/stable/visualization/matplotlib_integration.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⛏ Exercise 2.0\n",
    "\n",
    "- Create a `Table` called `tab1` with one column with values `[np.arange(10)]` named `[velocity]`.\n",
    "- Create a `QTable` called `qtab1` with one column with values `[np.arange(10)]` named `[velocity]`.\n",
    "- Show or print both tables. Do you see any difference?\n",
    "- Compute the mean  value of `tab1['velocity']` using `np.mean()` and also using `.mean()`. Do the same for `qtab1['velocity']`\n",
    "- Compute the square `()**2` of the velocity column, for both tables. What are the differences?\n",
    "- Execute `tab1['velocity'].data`. What do you obtain?\n",
    "- Do you obtain the same with `qtab1['velocity'].data`? Try `np.array(qtab1['velocity'])`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selection = tab['y'] > 20\n",
    "selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab[selection]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab['Date'][selection]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selection2 = (tab['x'] > 0) & (tab['MJD'] < 56400)\n",
    "selection2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab[['MJD', 'skycoord', 'x']][selection2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Storing tables on disk\n",
    "The Enhanced Character Separated Values table (`.ecsv`) format is a modified `csv` format that can handle metadata from astronomical tables and still keep the universaility and easy accessibility of `csv` files. It is able to store the data types, physical units and descriptions. It uses YAML, which is easy to read by humans and machines and can be opened with a plain text editor. It is not astronomy-specific. More details can be found in the [following link](https://github.com/astropy/astropy-APEs/blob/main/APE6.rst), including the comparison with other formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab.write('test_tab.csv', format='ascii.csv', overwrite=True)\n",
    "tab.write('test_tab.ecsv', format='ascii.ecsv', overwrite=True)\n",
    "tab.write('test_tab.html', format='jsviewer', overwrite=True)\n",
    "tab.write('test_tab.tex', format='ascii.latex', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All available formats are explained in the documentation `tab.write?`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load our problem data from last session"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⚠️  If you didn't create and save your own `data0` you can use the \"official\" data0 stored in the backup folder to be able to continue with this tutorial. Note that these commands will don't extract the \"official\" file if `data0.ecsv` already exists. If you want to discard your file and use the official one, first open a terminal and manually remove your `data0.ecsv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "if not os.path.isfile('../data/data0.ecsv'):\n",
    "    print('Using official data0 file')\n",
    "    os.system('unzip ../data/backup/data0.ecsv.zip')\n",
    "    os.system('mv data0.ecsv ../data')\n",
    "else:\n",
    "    print('Doing nothing because ../data/data0.ecsv already exists')           "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data0 = Table.read('../data/data0.ecsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data0.show_in_browser(jsviewer=True)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### ⛏ Exercise 2.1\n",
    "- Print the first 5 rows of the table with `data0[0:5]`.\n",
    "- Obtain a list of columns of the table `data0`.\n",
    "- Create a new table `my_table` containing only the columns Right Ascension, Declination, and the two proper motions, and the associated errors for all of them. You should obtain a table with 8 columns.\n",
    "- Convert it to a `QTable` using `my_qtable = QTable(my_table)`. Work with `my_qtable` for the rest of the exercise.\n",
    "- Obtain the description (column name, format, units and description) using `.info()` method.\n",
    "- Print the mean (`.mean()`) R.A. in degress and the median Declination in arcmin.\n",
    "- Print the standard deviation (`.std`) of the R.A. in arcseconds and of the Declination in arcmin.\n",
    "- Compute the minimum and maximum, `.min()`, `.max()`, proper motion in right ascension in units of mas/yr\n",
    "- Compute the minimum and maximum proper motion in declination in units of arcmin/day\n",
    "- Create a new column `pm_error` in `my_qtable` with the total uncertainty in the position that accounts for the quadratic sum (`np.sqrt(()**2 + ()**2)`) of the uncertainty of the proper motion components in units of arcmin per minute.\n",
    "- Create a filter to select stars with `pm_error` < 0.0005.\n",
    "  - How many stars comply that condition?\n",
    "  - What is the mean `ra` and mean `de` of that subgroup?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Additional fun\n",
    "\n",
    "Repeat all the previous steps but with `my_table` instead of `my_qtable`. What are the main limitations you find when working with `Table` units?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 2.2\n",
    "\n",
    "Compute a new column for `data0` with the absolute magnitude in the Gaia filter, g. Call the new column `Mg` and use the formula:  \n",
    "$M_{\\rm V} = m + 5log_{\\rm 10}p - 10$,  \n",
    "where $M_{\\rm V}$ is the absolute magnitude, $m$ is the apparent magnitude and $p$ is the stellar parallax in milliarcseconds.\n",
    "\n",
    "For example:\n",
    "\n",
    "Vega has a parallax $p$ of 0.129 arcsec, and an apparent magnitude $m_{\\rm V}$ of 0.03:  \n",
    "$M_{\\rm V} = 0.03 + 5log_{\\rm 10}129 - 10$ = 0.58\n",
    "\n",
    "The column containing the apparent magnitude in the Gaia band is called `phot_g_mean_mag`, and the parallax is `parallax`. The logarithm is computed as `np.log10(<relevant column>)`.\n",
    "\n",
    "You have more details, for instance, here: https://en.wikipedia.org/wiki/Absolute_magnitude#Examples\n",
    "\n",
    "- Print the first few rows of `data0` and make sure the new column is created.\n",
    "\n",
    "- Verify that your solution is right with this code:  \n",
    "`data0['Mg'][0] - 8.218870196400683 < 1e-5`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory analysis\n",
    "## Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "#from astropy.visualization import quantity_support   # If you use QTable, you will need this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(data0['ra'],     # X variable\n",
    "         data0['dec'],    # Y variable\n",
    "         '.',             # marker, can also be 'o', '+', '>', 's', etc..\n",
    "         color='k',       # color k:black, g:green, r:red, b:blue, etc\n",
    "         alpha = 0.2,     # transparency, between 0 and 1\n",
    "         ms=1);           # marker size "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not very informative. Let's create a `figure` and `axes` with `subplots`. Let's make the image bigger and use `scatter` to be able to set a symbol size `s` proportional to the flux of the star. We will also add some axes labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,10))\n",
    "\n",
    "# We use scatter here because it allows us to plot points with different sizes/colors\n",
    "ax.scatter(data0['ra'],\n",
    "           data0['dec'],\n",
    "           s = data0['phot_g_mean_flux']/1e5)  # Size of the data points, here proportional to the star flux\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "# Here we invert the direction of the right ascension axis\n",
    "ax.invert_xaxis()\n",
    "\n",
    "ax.set_xlabel('Right Ascension [deg]')  # Write X axis label\n",
    "ax.set_ylabel('Declination [deg]');     # Write Y axis label"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot any pair of variables. Here a quick plot without too many configuration for quick lookup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(data0['pmra'], data0['pmdec'], '.k')   # '.k' is used as an abbreviation for the format marker='.', color='k', ls='' (black points without lines)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we manually fix the X and Y limits so we can zoom in the relevant region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(data0['pmra'],\n",
    "        data0['pmdec'],\n",
    "        '.k')\n",
    "\n",
    "ax.set_xlim(-80, 60)\n",
    "ax.set_ylim(-80, 60);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also we can use log scale when that is relevant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(data0['phot_g_mean_flux'],\n",
    "        data0['parallax_error'], '.k', ms=1, alpha=0.5)  # ms is markersize and alpha is the transparency, between 0 and 1\n",
    "\n",
    "ax.loglog();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But! Never forget to add labels to identify what you are plotting!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(data0['phot_g_mean_flux'],\n",
    "        data0['parallax_error'], '.k', ms=1, alpha=0.5)\n",
    "\n",
    "ax.loglog()\n",
    "\n",
    "ax.set_xlabel('G-band mean flux [e/s]')\n",
    "ax.set_ylabel('Parallax error [mas]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tracking the axis label and units manually means that it is possible that we make mistakes. We can use the column information instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data0['phot_g_mean_flux'].description)\n",
    "print(data0['parallax_error'].unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⛏ Exercise 2.3\n",
    "- Repeat the plot above (phot_g_mean_flux vs parallax_error) but using the methods `.description` and `.unit` to write the x and y labels automatically.\n",
    "\n",
    "- Second, repeat the plot again but now make the script more generic: start by defining `col1` and `col2` with the names of the relevant columns and use only those variables so you don't need to write the explicit column name more than once:\n",
    "\n",
    "```python\n",
    "col1 = 'phot_g_mean_flux'\n",
    "col2 = 'parallax_error'\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⛏ Exercise 2.4\n",
    "- Repeat the plot but now use the columns `phot_g_mean_mag`, `parallax`.\n",
    "- Can you interpret the horizontal line with an overdensity of bright stars?\n",
    "- Compute the average parallax of the stars with magnitude < 10. Define a variable `selection`.\n",
    "- Repeat the plot but add one more `ax.plot()` to show in the X axis `data0[col1][selection]` and in the Y axis `data0[col2][selection]`. Make the symbols larger.\n",
    "- Use the `np.mean` and `np.median` functions to find the mean and median parallax of the selection.\n",
    "\n",
    "Tip: you may not be able to compute the **median** because there are missing values and it produces a masked array. You can use the method `.compressed()` to the column selection to be able to compute `np.median`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Additional fun\n",
    "Plot horizontal and vertical lines with ax.hline and ax.vline on the relevant magnitude and parallax computed in the previous exercise. You can use `ax.axvline` and/or `ax.axhline` functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coordinate plots in sky projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "ax = fig.add_subplot(111, projection=\"mollweide\")\n",
    "\n",
    "ax.scatter(data0['ra'].to(u.radian),\n",
    "           data0['dec'].to(u.radian), marker='.', color='r')\n",
    "\n",
    "ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's overplot the galactic plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galactic_plane = SkyCoord(l=np.arange(-180, 180),\n",
    "                          b=np.zeros(360),\n",
    "                          frame='galactic', unit=u.deg)\n",
    "\n",
    "galactic_plane_eq = galactic_plane.transform_to('icrs')\n",
    "gal_ra  = galactic_plane_eq.ra.wrap_at('180d').radian\n",
    "gal_dec = galactic_plane_eq.dec.radian\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "ax = fig.add_subplot(111, projection=\"mollweide\")\n",
    "ax.scatter(data0['ra'].to(u.radian), data0['dec'].to(u.radian), marker='.', color='r')\n",
    "plt.plot(gal_ra, gal_dec, 'k.')\n",
    "\n",
    "ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Additional fun\n",
    "Make the previous plot but in Galactic coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Additional fun\n",
    "Prepare a plot `pmra` vs `pmdec` including errorbars using the matplotlib function `plt.errorbars`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histogram distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=1, figsize=(10,8))\n",
    "\n",
    "ax.hist(data0['parallax'], bins=np.arange(-5, 15, 0.5))\n",
    "\n",
    "ax.set_xlabel('Parallax [mas]')\n",
    "ax.set_ylabel('Number of stars');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### ⛏ Exercise 2.5\n",
    "\n",
    "Modify the line `ax.set_xlabel('Parallax [mas]')` to automatically find the units of the column being plotted. Do the same we did before when we defined the function `str_label()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the same histogram as before but starting with the generic variable col1. Do not use 'parallax' in any other part\n",
    "col1 = 'parallax'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 2.6\n",
    "\n",
    "First of all, we see that there are negative parallaxes, which don't have physical meaning, but are a consequence of errors in the parallax determination. We can create a second `Table` named `data1` that ignores any negative parallax. We will use this table from now on.\n",
    "\n",
    "- Define the variable `positive_parallaxes` as those entries with parallax greater than 0. `data0['parallax'] > 0`. This is a boolean array with `True` values for rows with positive parallaxes, and `False` otherwise.\n",
    "- Now slice table `data0` by selecting the rows with `True` values: `data0[positive_parallaxes]`. Assign that new table to `data1`.\n",
    "- Print the lenght of `data0` and `data1`. How many stars each one has? Use `len(data0)`, `len(data1)`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equivalencies\n",
    "Let's work with distances in kpc, which are more familiar to us than parallaxes. We will use an astropy unit transformation as before. However, an angle (mas) cannot be converted to distance (kpc) without knowing how the transformation should occur. We need to parse which equivalency to use to make the transormation.\n",
    "\n",
    "There is a lot of information in section [Equivalencies](https://docs.astropy.org/en/stable/units/equivalencies.html), for example to convert spectral units (nm to Hz) or conversions from wavelength/frequency/energy including doppler effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1['distance'] = data1['parallax'].to(u.kpc, equivalencies=u.parallax())\n",
    "data1['distance'].description = 'Distance from Earth'\n",
    "data1[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### ⛏ Exercise 2.7\n",
    "\n",
    "Compute manually the distance in kpc and check that the transformation has worked. Start by using the numpy array without units `np.array(data1['parallax'])` and compute the distance in kpc  with the formula $d [{\\rm kpc}]= \\frac{1}{p[{\\rm mas}]}$, where $p$ is the parallax. Compute the average of the residual `data1['distance'] - d_kpc`, where `d_kpc` is the manually computed distance in units of `kpc`. What is the average residual?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 2.8\n",
    "- Going back to the parallax distribution. Plot again the histogram of the `data[parallax]`, starting with `col1='parallax'`. Make the figure size slightly bigger, and use a smaller bin size. Try different `bins` selection either with `np.arange` (which fixes the step) or `np.linspace` (which fixes the number of steps) until you see some significant structure. Identify, approximately, the range of parallaxes that look interesting to you. (You can use `ax.set_xlim()` to tune the plot range).\n",
    "- Make a similar plot but with `col1='distance'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see a very interesting accummulation of stars at a parallax of approximately 5.2 mas. We can create a filter to select the start in that particular range. We will overplot the distribution of the whole sample and the one of the selected group. Check how many stars we have selected with the filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "manual_filter1 = (data1['parallax'] > 5.0*u.mas) & (data1['parallax'] < 5.7*u.mas)\n",
    "\n",
    "cluster1 = data1[manual_filter1]\n",
    "cluster1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an alternative way to deal with this problem. Instead of creating a new `Table` we could create a new Boolean column with True/False that indicates if each star is part of the cluster or not. We will not use this new column `cluster`, but it is good so we have all the information together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1['cluster'] = manual_filter1\n",
    "data1[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The get the stars from the cluster we can simply select the rows matching the filter, i.e., selecting the `True` values. We could use `data1[manual_filter1]`, but we do not want to carry the filter around, it is better if the information is embedded in the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1[data1['cluster']]  # This is exactly the same as cluster1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will not use the 'cluster' column. It is created here just for training purposes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 2.9\n",
    "- Use the last two plots with the parallax and the distance distributions. Appart from the `data1`, add now the histograms for the table `cluster1`. You should see the two distributions overlapped.\n",
    "- Set the labels to 'Full sample' and 'Cluster', respectively.\n",
    "- Fine tune the data ranges as needed.\n",
    "- What is the approximate distance of the selected cluster?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 2.10\n",
    "\n",
    "Create the `phot_g_mean_mag` vs `parallax_error` we did before, but including data from the two tables: `data1` in black and `cluster1` in red."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the cluster is dominated by bright stars (lower magnitude), and form a sharp cluster in parallax/distance, as we already knew."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial distribution of the cluster\n",
    "We plot the distribution of start in the sky. First, all the stars in the sample are plotted in grey. The stars of the cluster are plotted in color, with the colorscale representing the distance from the Earth in pc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,10))\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.scatter(data1['ra'], data1['dec'], c='gray', s=1, alpha=0.5);\n",
    "l = ax.scatter(cluster1['ra'], cluster1['dec'], c=cluster1['distance']*1000., s=20);  # We save the variable l to be used in the colorbar\n",
    "\n",
    "ax.set_xlabel('Right Ascension [deg]')\n",
    "ax.set_ylabel('Declination [deg]');\n",
    "\n",
    "# Here we invert the direction of the right ascension axis\n",
    "ax.invert_xaxis()\n",
    "\n",
    "# Show the color bar\n",
    "cb = fig.colorbar(l);\n",
    "cb.set_label('Distance [kpc]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's include the photometry information and use a different color for the `cluster1` stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=1, figsize=(15,10))\n",
    "ax.scatter(data0['ra'], data0['dec'], s=data0['phot_g_mean_flux']/1e5);\n",
    "ax.scatter(cluster1['ra'], cluster1['dec'], s=cluster1['phot_g_mean_flux']/1e5);\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "# Here we invert the direction of the right ascension axis\n",
    "ax.invert_xaxis()\n",
    "\n",
    "ax.set_xlabel('Right Ascension [deg]')\n",
    "ax.set_ylabel('Declination [deg]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 🌪 Additional fun\n",
    "Repeat the previous plot by making it more general, using generic columns `col1`, `col2`, and automatic labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no apparent pattern of the selected stars, although there seems to be an overdensity at the center, specially in Right Ascension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ra_range = [np.min(data1['ra']), np.max(data1['ra'])]\n",
    "                 \n",
    "                                       \n",
    "ax.hist(data1['ra'],    bins=np.arange(ra_range[0], ra_range[1], 0.04), label='Full sample')\n",
    "ax.hist(cluster1['ra'], bins=np.arange(ra_range[0], ra_range[1], 0.12), label='Cluster')\n",
    "\n",
    "ax.set_xlabel('Right Ascension [deg]')\n",
    "ax.set_ylabel('Number of stars');\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots( figsize=(8,6))\n",
    "\n",
    "de_range = [np.min(data1['dec']), np.max(data1['dec'])]\n",
    "\n",
    "ax.hist(data1['dec'],    bins=np.arange(de_range[0], de_range[1], 0.02), label='Full sample')\n",
    "ax.hist(cluster1['dec'], bins=np.arange(de_range[0], de_range[1], 0.08), label='Cluster')\n",
    "\n",
    "ax.set_xlabel('Declination [deg]')\n",
    "ax.set_ylabel('Number of stars');\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nothing obvious or much interesting. The overdensity in the centre for R.A. is clear, but in Dec it is not clear, there may be a slope."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What it is clear is that the cluster is not easily identified in R.A. and Declination. We will need to explore more variables to characterize the cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✨ Exercise 2.11\n",
    "\n",
    "Let's save the data sets for future use. Like in session 1, we will save them with `format='ascii.ecsv'`\n",
    "- Use method `data1.write()` to save the table in file `'../data/data1.ecsv'` with the correct format.\n",
    "- Use method `cluster1.write()` to save the table in file `'../data/cluster1.ecsv'` with the correct format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pandas tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = data1.to_pandas()\n",
    "df1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring a pandas table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first thing we see is that we have lost the unit information. That is an important problem if we don't track the column operations properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1['ra']\n",
    "df1[['ra','dec']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1['pmra'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To slice by index value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.loc[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To slice by index number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.iloc[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can be combined with slices in one or several columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1['pmra'].iloc[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.iloc[5:8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1[['pmra','pmdec']].iloc[6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1[['pmra','pmdec']].iloc[6:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pandas operations. Aggregate and groupby"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily do operations to columns. First by aggregating values according to some functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.aggregate(['sum', 'min'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.aggregate({'ra' : ['mean', 'min', 'max', 'std'],\n",
    "               'dec' : ['mean', 'min', 'max', 'std'],\n",
    "               'parallax': 'std'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_func(row):\n",
    "    if row['cluster']:\n",
    "        error =  np.sqrt(row['pmra_error']**2 + row['pmdec_error']**2)\n",
    "    else:\n",
    "        error = None\n",
    "    return error\n",
    "\n",
    "df1['pm_error'] = df1.apply(lambda row: my_func(row), axis=1)\n",
    "\n",
    "df1[df1['cluster']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1[~df1['cluster']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NaN values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1['dr2_radial_velocity'].isna().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_rv = df1[['dr2_radial_velocity', 'parallax']].dropna()\n",
    "df_rv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Groupby"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.groupby('cluster').aggregate('mean')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster1['astrometric_matched_transits'].description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1['astrometric_matched_transits'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make many groups according to how many observations a star has."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.groupby('astrometric_matched_transits').aggregate(['mean','min'])[['parallax_error','pmra_error', 'pmdec_error']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.groupby('astrometric_matched_transits').aggregate(['mean','min'])[['parallax_error']].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working with datetime series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "times = pd.date_range(\"2021-06-9\", periods=N, freq='D')\n",
    "\n",
    "ts = pd.DataFrame({'v1': np.random.normal(0.5,100,N) + 0.1*np.arange(N)**2,\n",
    "                   'v2': np.random.normal(0.5,100,N) - 0.1*np.arange(N)**2},\n",
    "                   index=times)\n",
    "ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.plot.scatter(x='phot_g_mean_flux', y='parallax_error',  marker='.', s=0.5, alpha=0.1)\n",
    "plt.loglog();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(data0['phot_g_mean_flux'], data0['parallax_error'], '.k', ms=1, alpha=0.1)\n",
    "plt.loglog();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1[['pmra','pmdec']].plot.hist(bins=np.arange(-80, 60, 3), alpha=0.5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1[['phot_g_mean_mag','phot_bp_mean_mag','phot_rp_mean_mag']].plot.box();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1[['phot_rp_mean_mag', 'phot_bp_mean_mag', 'phot_g_mean_mag']].plot.hist(bins = 80, alpha=0.8);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.plot.hexbin(x='ra', y='dec', gridsize=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas.plotting import scatter_matrix\n",
    "scatter_matrix(df1[['pmra_error','pmdec_error','pmra','pmdec']], alpha=0.2, diagonal='kde', figsize=(8,8));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many more examples here: https://pandas.pydata.org/docs/user_guide/visualization.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}