pandas_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Everyone should use Pandas! (and sometimes astropy tables)\n",
    "* This tutorial is meant to get you started with using Pandas and Astropy tables to organize and work with your data sets \n",
    "\n",
    "* Storing your data as pandas dataframes makes querying and sorting your data, merging different data sets, doing statistics, and visualiziing your data faster, easier, and cleaner. \n",
    "* Pandas also handels messy data and data with missing fields very well \n",
    "* Pandas is a python package widely used to build and work with dataframes. Astropy has a class called Tables which has similar functionality but much more limited and less supported. I introduce it breifly here because it has a few astronomy specific features absent in Pandas that are sometimes useful. For most things the two objects can be used interchangably. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Main documentation for both Pandas, Astropy tables, and Astroquery\n",
    "* [Pandas](http://pandas.pydata.org/pandas-docs/stable/)\n",
    "* [Astropy Tables](http://docs.astropy.org/en/stable/table/)\n",
    "* [Astroquery](https://astroquery.readthedocs.io/en/latest/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#import all needed python packages\n",
    "from __future__ import print_function\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "#comment out the next two lines if you do not have astroquery installed\n",
    "from astroquery.sdss import SDSS\n",
    "from astroquery.vizier import Vizier\n",
    "\n",
    "from astropy.table import Table, Column, join\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Types of Tables in Pandas\n",
    "1. Series - 1D array, can be build out of a python dictionary \n",
    "2. Dataframes - 2D array "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Tables are different from numpy arrays:\n",
    "* Each row in a table is assigned an index. \n",
    "* In a dataframe column also has a name \n",
    "* A numpy array can also be indexed but you can't change what the indecies are named \n",
    "    (always indexed in the same way)\n",
    "* tables can have columns of all different data types "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np_array = np.array([1,2,3,4,5])\n",
    "\n",
    "print('Numpy Array')\n",
    "print(np_array)\n",
    "print('\\n')\n",
    "\n",
    "#It will assign the index for you if you don't give it one\n",
    "#Indexing always starts at 0 in python\n",
    "pd_series = pd.Series([1,2,3,4,5])\n",
    "\n",
    "print('Pandas Series')\n",
    "print(pd_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Series in Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#You can define the index for the series \n",
    "s1 = pd.Series([1,2,3,4,5],index = ['a','b','c','d','e'])\n",
    "\n",
    "print(s1, '\\n')\n",
    "\n",
    "#building a series from a dictionary \n",
    "data = {'cats' : 15., 'dogs' : 35., 'turtles' : 100.}\n",
    "s2 = pd.Series(data)\n",
    "\n",
    "print(s2, '\\n')\n",
    "\n",
    "#You can assign the index (pandas handels missing data very well)\n",
    "s3 = pd.Series(data,index=['turtles','cats','birds','dogs'])\n",
    "\n",
    "print(s3, '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Dataframes in Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#example of how to build a dataframe\n",
    "#Example: I have an array of data that tells me the number of animals on a given date\n",
    "dates = pd.date_range('20180620', periods=6)\n",
    "rand_data = np.random.randint(low=0, high=100, size=(6,4))\n",
    "df = pd.DataFrame(rand_data, index=dates, columns=['cats','dogs','birds','turtles'])\n",
    "\n",
    "print(df)\n",
    "print('\\n')\n",
    "\n",
    "#If you don't care what the index is you don't have to give it one\n",
    "#If you don't assign and index it just gives them index values from 0 to # rows\n",
    "\n",
    "df = pd.DataFrame(rand_data, columns=['cats','dogs','birds','turtles'])\n",
    "\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Importing data with Pandas \n",
    "* Pandas can easily read in files of data in one line \n",
    "* Unless you specify the dtypes you want for each column it chooses the best match \n",
    "* Each read function in Pandas is very flexible with a lot of arguments to to read in a lot of differnet formats\n",
    "* Generic delimited text [(documentation)](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_table.html)\n",
    "    * `read_csv()` \n",
    "    * `read_table()` - defaults to tab separated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Other formats\n",
    "Pandas has built-in support for [reading from or to other common formats/sources](http://pandas.pydata.org/pandas-docs/stable/io.html):\n",
    "* Excel -- `read_excel()`\n",
    "* JSON -- `read_json()`\n",
    "* SQL -- `read_sql()`\n",
    "* Stata -- `read_stata()`\n",
    "* SAS (XPORT or SAS7BDAT) -- `read_sas()`\n",
    "* etc..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## We are going to work with the HETDEX Pilot Survey emission line catalog\n",
    "* this is a Vizier catalog that can be found [here](http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/ApJS/192/5/hetdex)\n",
    "* The page has an option to save the table as a tab separated file which I did "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#importing data with pandas \n",
    "\n",
    "#Because the file has such a large header you need to tell it which row you want it to use for the column names (header)\n",
    "#Tell it to skip the rows below it that have the units and the line \n",
    "#If you don't care about the header information it is sometimes easier to just delete it from the top of the file\n",
    "#  if you delete the header just make sure to leave the top row of column names\n",
    "hps_pd = pd.read_table('hps_catalog.txt', delimiter=None, header=62, skiprows=[64,65])\n",
    "hps_pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Astropy supports different file formats (commonly used in astronomy) that Pandas doesn't\n",
    "* some of the formats are more astronomy specific (ex. fits)\n",
    "* It does allow you to save tables as deluxe latex tables which is extremely useful!\n",
    "* It also will read in a row of units for columns (useful when reading in catalogs)\n",
    "* A [list of the format Astropy supports](http://docs.astropy.org/en/stable/io/unified.html#built-in-readers-writers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#importing data with astropy tables\n",
    "\n",
    "hps_at = Table.read('hps_catalog.txt', format='ascii')\n",
    "hps_at\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Reading in with Pandas vs. Astropy Tables \n",
    "* You can see that Tables did not automatically pick up on there being two rows under the column names that are not data\n",
    "* the Tables read function does not have arguments for formating files like pandas \n",
    "* so unless your data fits the formats provided astropy tables can be annoying to work with \n",
    "* I usually have to format the file before I try to read it in with astropy tables \n",
    "* Pandas provides a lot more flexlibilty "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## However, you can also load it in with Astropy Tables using Astroquery\n",
    "* Astroquery allows you to load in a catalog from a bunch of different sources \n",
    "* [Here](https://astroquery.readthedocs.io/en/latest/#catalog-archive-and-other) is a list of all of the catalogs and archives that now work with astroquery (some are better supported than others)\n",
    "* [Here](http://astroquery.readthedocs.io/en/latest/gallery.html) is a page with astroquery examples for some of the databases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#queries with astropy (Vizier)\n",
    "\n",
    "hps_catalog_list = Vizier.find_catalogs('J/ApJS/192/5/')\n",
    "hps_catalogs = Vizier.get_catalogs(hps_catalog_list.keys())\n",
    "hps_at = hps_catalogs[0]\n",
    "\n",
    "hps_at\n",
    "#hps_at.show_in_notebook() #This is a tool for viewing astropy tables in a Jupyter notebook\n",
    "#hps_at.show_in_browser(jsviewer=True) #this is a tool for viewing a searchable astropy table in a web browser"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Interfacing between pandas and astropy tables \n",
    "* you can covert from one to the other very easily with the following commands "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# astropy tables to pandas\n",
    "hps_pd_from_at = hps_at.to_pandas()\n",
    "\n",
    "# pandas to astropy table\n",
    "hps_at_from_pd = Table.from_pandas(hps_pd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# **From now on in this tutorial we will be using Pandas**\n",
    "* some of this functionality exists in Astropy Tables but not all of it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exploring your data set "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#hps_pd.head()\n",
    "hps_pd.tail(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "hps_pd.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "hps_pd.columns.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "hps_pd.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "hps_pd.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Indexing Pandas dataframes\n",
    "* You can either index by the row and column number as you with in numpy \n",
    "* Or more commonly you use the columns names and the row index value\n",
    "* Pandas defults to using label indexing\n",
    "* For clarity pandas makes you specify if you using integer or label indexing\n",
    "    * this is neccessary because the columns names are allowed to be numbers \n",
    "    * without specifying pandas would not know if you mean column number 1 or the column named 1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#label indexing\n",
    "\n",
    "#hps_pd['EW']\n",
    "#hps_pd.EW\n",
    "hps_pd.loc[:,'EW']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#You can ask for certian rows in that column\n",
    "hps_pd['EW'][2:5] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#can also get more than one column\n",
    "hps_pd.loc[:25,['Trans','Lya','S/N']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#I want all of the detections at a z<0.5 with and EW>50\n",
    "\n",
    "#You can also use the query function\n",
    "#hps_pd.query('(z < 0.5) & (EW > 50.)')\n",
    "\n",
    "#can also write it this way \n",
    "hps_pd.loc[(hps_pd.z < 0.5) & (hps_pd.EW > 50)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#You can also index by location like you do in numpy using iloc\n",
    "\n",
    "#If you give it a single value it returns a series of that row number \n",
    "hps_pd.iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# You can also give it a specific section of the table you want \n",
    "#format is [row_start:row_end, col_start:col_end]\n",
    "hps_pd.iloc[0:10, 2:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#It will also take a list of column numbers \n",
    "hps_pd.iloc[0:10, [1,4,6,9]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Useful Functions\n",
    "* Pandas has a lot of useful built in functions \n",
    "* Pandas apply function allows you to operate a function on all values in a column/table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "hps_pd['Trans'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#You may notice when it read in the Trans column it added extra spaces making it difficult to query\n",
    "# You can use the apply function to strip each value in the Trans column \n",
    "# Anonymous (lambda) funcitons in python are a convienent way to write a one line function to operate over a variable\n",
    "# the strip() function strips a string of its whitespace\n",
    "hps_pd['Trans'] = hps_pd['Trans'].apply(lambda x: x.strip())\n",
    "\n",
    "#can ues the values function to return an array of just a columns values \n",
    "hps_pd['Trans'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#find the unique types of emission line types \n",
    "emission_lines = hps_pd['Trans'].unique()\n",
    "emission_lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#get an table of HPS ids, fluxes, and S/N for Ly{alpha} lines where the probablilty of it being Lya is greater than 50%\n",
    "hps_pd.query('(Trans == \"Ly{alpha}\") & (Lya > 0.5)').loc[:,['HPS','Flux','S/N']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# you want the lya line with the largest flux \n",
    "max_lya_flux = hps_pd.query('(Trans == \"Ly{alpha}\") & (Lya > 0.5)')['Flux'].max()\n",
    "print('MAX Flux: ', max_lya_flux, '\\n')\n",
    "\n",
    "max_flux_obj = hps_pd.loc[hps_pd['Flux'] == max_lya_flux]\n",
    "print(max_flux_obj, '\\n')\n",
    "#hps_pd.iloc[max_lya_flux]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#want to know how many of each line\n",
    "hps_pd['Trans'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#count the number of LyA lines with a probablility greater than 50%\n",
    "hps_pd.query('(Trans == \"Ly{alpha}\") & (Lya > 0.5)')['HPS'].count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#sort the LyA lines by the flux values in descending order \n",
    "LyA_sorted_by_flux = hps_pd.query('Trans == \"Ly{alpha}\"').sort_values('Flux', ascending=False)\n",
    "LyA_sorted_by_flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#build function to give to apply for a column\n",
    "def AA_to_nm(wl):\n",
    "    wl = wl/10\n",
    "    return wl\n",
    "\n",
    "hps_pd['lambda'].apply(AA_to_nm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Here is a function that builds all of the RAj2000 and DEJ2000 values into actual astropy coordinates\n",
    "#It stores the coordinates in a new column in the dataframe called coord\n",
    "\n",
    "def build_astropy_coords(row):\n",
    "    ra  = row['RAJ2000'].split()\n",
    "    dec = row['DEJ2000'].split()\n",
    "    c = SkyCoord(ra[0]+'h'+ra[1]+'m'+ra[2]+'s', dec[0]+'d'+dec[1]+'m'+dec[2]+'s', frame='icrs')\n",
    "    return c\n",
    "\n",
    "#You give apply the function and the argument axis=1 to it feeds the function the rows of the dataframe\n",
    "#It stores these new values in a new column called coord\n",
    "hps_pd['coord'] = hps_pd.apply(build_astropy_coords, axis=1)\n",
    "\n",
    "#prints the new ra and dec of the stored coordinate for the first row\n",
    "#print(hps_pd['coord'][0].ra, hps_pd['coord'][0].dec)\n",
    "\n",
    "hps_pd.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#want to find all of the rows that have an x-ray detection \n",
    "\n",
    "#first want to replace all of the rows with nothing but spaces in the X-ray column with nan values \n",
    "hps_pd['X-ray'] = hps_pd['X-ray'].replace('                ', np.nan, regex=True)\n",
    "\n",
    "#use the drapna function to remove the rows without a detection. \n",
    "#the .index function returns the indecies of the columns with an x-ray detection\n",
    "#Then use the .loc function to find the rows with that index \n",
    "possible_AGN = hps_pd.loc[hps_pd['X-ray'].dropna().index]\n",
    "possible_AGN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#the dropna function had many arguments that make it very powerful. \n",
    "\n",
    "hps_drop_allNaN  = hps_pd.dropna()     #drop all rows that have any NaN values\n",
    "hps_drop_onlyNaN = hps_pd.dropna(how='all')     #drop only if ALL columns are NaN\n",
    "hps_drop_someNaN = hps_pd.dropna(thresh=3)   #Will only keep rows with 3 or more **not** NaN values\n",
    "\n",
    "#this is a better way of writing the line defining possible_AGN in the previous cell\n",
    "possible_AGN = hps_pd.dropna(subset=['X-ray'])   #Drop only if NaN in specific column \n",
    "possible_AGN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#changing/appending data \n",
    "#saw already how to rename columns but here is how to add/remove columns "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Building groups\n",
    "* you can build groups based on column values \n",
    "* This could be really useful in this example to group by the line type (Trans column)\n",
    "* pandas provides many built-in operations on groups--e.g., mean, count, etc.\n",
    "* see http://pandas.pydata.org/pandas-docs/stable/groupby.html#groupby-object-attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Split the dataset on Trans. This creates a pandas GroupBy object.\n",
    "lines = hps_pd.groupby('Trans')\n",
    "\n",
    "# You can call one of the groups by the name\n",
    "NeV_lines = lines.get_group(\"[NeV]3426\")\n",
    "print(NeV_lines)\n",
    "\n",
    "lines['Flux'].mean()\n",
    "\n",
    "# The above line can also be more explicitly written as:\n",
    "# groups['age'].aggregate('mean')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Types of joins\n",
    "* Joins allow you to join together different dataframes based on a key values (similar to SQL)\n",
    "* for more information see the [merging documentation](https://pandas.pydata.org/pandas-docs/stable/merging.html)\n",
    "* `left`  joins - Use keys from left frame only\n",
    "* `right` joins - Use keys from right frame only\n",
    "* `outer` joins - Use union of keys from both frames\n",
    "* `inner` joins - Use intersection of keys from both frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Diagram of joins\n",
    "<img src=\"https://i.stack.imgur.com/hMKKt.jpg\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Load in new a table from a paper to merge with the HPS catalog\n",
    "#Bridge et. al. 2014 derived properties from SED fitting for the OII emitters in the HPS catalog\n",
    "\n",
    "bridge_cols=['HPS_name', 'logM', 'logSFR', 'E(B-V)', 'OII_z']\n",
    "bridge = pd.read_table('bridge2015_cat.ascii', delim_whitespace=True, names=bridge_cols, skiprows=16)\n",
    "bridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#We want to join this table with the HPS catalog\n",
    "#They share the common key (HPS id number) but they have different names in each table\n",
    "\n",
    "#We first change the name of the HPS id in the bridge dataframe\n",
    "bridge.rename(columns={'HPS_name':'HPS'}, inplace=True)\n",
    "\n",
    "#We can then merge the two dataframes on the key HPS\n",
    "#We use and inner join so that we only get rows back that have both HPS and bridge data\n",
    "#If we use any other join it would fill in the missing data with NaNs \n",
    "HPS_OII = pd.merge(hps_pd, bridge, how='inner', on=['HPS'])\n",
    "HPS_OII"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Here is an example of a outer join\n",
    "* In this case we want to preserve all of the rows in the HPS catalog \n",
    "* but when we can we want to add a value for the field that emission line was found in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read in the table from the vizier catalog that contains the information on the fields observed \n",
    "hps_fields = pd.read_table('hps_fields.txt', delimiter=None, header=50, skiprows=[52,53])\n",
    "print(hps_fields)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#this function just finds the hour of the RA in the column 'RAJ2000'\n",
    "def find_RA_hour(row):\n",
    "    ra_h = row['RAJ2000'].split()[0] \n",
    "    return ra_h\n",
    "\n",
    "#apply this function to build an RA_hour column in the hps dataframe\n",
    "hps_pd['RA_hour'] = hps_pd.apply(find_RA_hour, axis=1)\n",
    "print(hps_pd['RA_hour'].unique())\n",
    "\n",
    "#apply this function to build an RA_hour column in the hps_fields dataframe\n",
    "hps_fields['RA_hour'] = hps_fields.apply(find_RA_hour, axis=1)\n",
    "\n",
    "#we only want to add the field name to the hps dataframe so this is just the Field column and the RA_hour to merge on\n",
    "fields = hps_fields.loc[:,['Field','RA_hour']]\n",
    "\n",
    "#We can then join the two databases to add the field column to the hps dataframe\n",
    "hps_join_fields = pd.merge(hps_pd, fields, how='outer', on=['RA_hour'])\n",
    "hps_join_fields.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Creating crosstabs \n",
    "* This is a useful tool that lets you analzye columns against eachother "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#You want to know how many of each line is in each field\n",
    "\n",
    "#creating crosstabs from groups\n",
    "ctab = pd.crosstab(hps_join_fields['Trans'], hps_join_fields['Field'])\n",
    "\n",
    "#You can also normalize the crosstab\n",
    "#ctab = pd.crosstab(hps_join_fields['Trans'], hps_join_fields['Field'], normalize='index')\n",
    "\n",
    "order = ctab.sum(1).sort_values(ascending=False).index\n",
    "ctab = ctab.loc[order, :]\n",
    "\n",
    "ctab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Saving dataframes \n",
    "* This can be done in one line with both pandas and astropy tables \n",
    "* For the most part the format of the output does not matter since it is just going to be loaded in as a dataframe later\n",
    "* in Pandas you can save the dataframe as a .csv or pickle file (may be more efficient for large files) \n",
    "* There are a lot more options for astropy tables \n",
    "    * basically any file you can load in you can save that file as the any of those outputs\n",
    "    * the most useful thing is you can save a datafame as a latex deluxe table "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#saving dataframes in Pandas \n",
    "HPS_OII.to_csv('bridge_HPS_cat.csv')\n",
    "\n",
    "#saving astropy table as a deluxe latex table \n",
    "#Here is an example of how we could make a table of just the possible AGN\n",
    "hps_at_OII = hps_at[hps_at['Trans'] == '[OII]']\n",
    "hps_at_OII.write('HPS_OII_table.tex', format='ascii.aastex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Pandas makes it easy to do statistics on columns and the entire dataframe\n",
    "* The describe function we saw earlier is customizable "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#statistics with panadas \n",
    "\n",
    "print(hps_pd.query('Trans == \"[OII]\"')['EW'].mean())\n",
    "print(hps_pd.query('Trans == \"[OII]\"')['EW'].std())\n",
    "\n",
    "# Describe the continuous variables in the dataset.\n",
    "# We can also include categoricals, or select other percentiles.\n",
    "#hps_pd.describe(percentiles=[0.1, 0.9]).round(2)\n",
    "\n",
    "#lines.get_group('[OII]').describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Data visualization using Pandas dataframes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Easy plotting with Matplotlib using Pandas columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#It is easy to plot in matplotlib by giving it pandas columns \n",
    "\n",
    "# Scatterplot of mass vs. SFR for the OII emitters from the Bridge et. al. 2014 paper\n",
    "plt.scatter(HPS_OII['logM'], HPS_OII['logSFR'], color='green', s=30)\n",
    "plt.xlabel(\"log(M)\")\n",
    "plt.ylabel(\"log(SFR)\")\n",
    "plt.title(\"Mass vs. SFR for the OII emitters in HPS\", fontsize=16);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plotting in Seaborn with dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#The crosstab we made can be easily plotted as a headmap \n",
    "ax = sns.heatmap(ctab[:12], annot=True, fmt='d', linewidths=.5,\n",
    "                 annot_kws={'size': 14}, cmap='Blues')\n",
    "\n",
    "# Customize the label size and figure size in base matplotlib.\n",
    "# Notice we can access all axis and figure properties through\n",
    "# # the handle we saved when calling seaborn.\n",
    "ax.tick_params(labelsize=14)\n",
    "ax.figure.set_size_inches((7, 8))  # update the plot size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot a heatmap of pairwise correlations between all continuous variables\n",
    "r = HPS_OII.corr().round(2)\n",
    "\n",
    "with sns.plotting_context(\"notebook\", font_scale=1.1):\n",
    "    ax = sns.heatmap(r, annot=True)\n",
    "    ax.figure.set_size_inches((12, 9.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#pairplots from a select group of columns for the OII catalog\n",
    "\n",
    "tmp_data = HPS_OII.loc[:,['Flux', 'logM', 'logSFR', 'E(B-V)','z']]\n",
    "\n",
    "with sns.plotting_context(\"notebook\", font_scale=1.3):  # Make text bigger\n",
    "    sns.pairplot(tmp_data, diag_kind='hist', kind='reg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.jointplot(x=\"logM\", y=\"logSFR\", data=HPS_OII, kind=\"kde\", color=\"purple\")\n",
    "g.plot_joint(plt.scatter, c=\"white\", s=30, linewidth=1, marker=\"+\")\n",
    "g.ax_joint.collections[0].set_alpha(0)\n",
    "g.set_axis_labels(\"$logM$\", \"$logSFR$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(6, 6))\n",
    "cmap = sns.cubehelix_palette(as_cmap=True, dark=0, light=1, reverse=True)\n",
    "g = sns.kdeplot(HPS_OII.logM, HPS_OII.logSFR, cmap=cmap, n_levels=60, shade=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plotting in Pandas "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HPS_OII.plot.scatter(x='logM', y='logSFR')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HPS_OII.plot.scatter(x='logM', y='logSFR', c='E(B-V)', s=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#Only want to plot histogram for lines with more than 15 detections\n",
    "many_detections = hps_pd.groupby('Trans').filter(lambda x: len(x)>15)\n",
    "\n",
    "# Plot flux histograms for all line types--with shared x-axis\n",
    "many_detections['Flux'].hist(by=many_detections['Trans'], bins=20, sharex=True);\n",
    "\n",
    "# Embiggen the figure a bit\n",
    "plt.gcf().set_size_inches((10, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_counts = hps_pd['Trans'].value_counts()\n",
    "print(line_counts)\n",
    "\n",
    "many_.plot.pie(subplots=True, autopct='%.2f', fontsize=15, figsize=(8, 8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}