Visualizing GAMA.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:01dfe52ee4e47fc834ae2f24c3866dd9898d6056120e692ec60f65fcc59471b3"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Acquiring and visualizing an extragalactic dataset in WorldWide Telescope\n",
      "using python, astroquery, astropy and pyWWT"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<img src=\"https://raw.githubusercontent.com/IPSScienceVisualization/python-tutorials/master/images/WWTbanner.PNG\">"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this tutorial we'll step through the entire process of pulling the data fom the GAMA survey database, performing some astrophysical calculations and sending the data to WWT."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Python"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Python is rapidly becomming the defacto astronomical programming language, supported by an active community of tool builders. Python's capabilities are extended through packages, and there are a number of good ones for dealing with astronomical data. This tutorial uses astropy which contains core astronomical functions and astroquery which lets you access astronomical databases from within the python environment. We also give an example of using the pyWWT package to directly interact with Worldwide Telescope.\n",
      "\n",
      "One of the challenges of working in python has always been installing all the packages needed to do your work. That has gotten much easier recently through some excellent package managers. If you are just starting out with python I strongly recommend using the Anaconda distribution from Continuum Analytics. Anaconda installs most of the packages needed for scientific computing with python (astropy is included with Anaconda, unfortunately you have to install astroquery and pyWWT yourself).\n",
      "\n",
      "This tutorial is provided as iPython notebook. You are either reading the notebook or an html rendering of the notebook. If the latter is the case you can download the notebook itself by clicking on the icon in the upper right corner of the webpage."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Let's get started - Step #1 Acquiring the GAMA Data from the database"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to get our galaxy catalog from the GAMA survey, to do that we need to query its database. A survey database typically contains several different tables each measuring a number of different quantities. Looking at the structure of the GAMA database (its schema) I see that the quantities I need to locate the galaxies in 3D space can be found in a table called specObj. These are sky position (RA and dec) and the redshift(z) which we will use to calculate the distance to the galaxy. There is another table called StellarMasses which has several interesting parameters, we will get three that might be useful to visualize later. There are the mass of the galaxy (logmstar), its instrinsic brightness (absMag_r) and its metalicity (metal).\n",
      "\n",
      "The database can be accessed using the Structured Query Language(SQL). SQL has a SELECT/FROM/WHERE format, after SELECT you list what parameters you want to pull from the dataase, after WHERE you specifiy which database tables they are located in, and after WHERE you can place conditions on which objects you wish to get back.\n",
      "\n",
      "So our query will be:\n",
      "\n",
      "SELECT SpecObj.ra, SpecObj.dec, SpecObj.z,StellarMasses.logmstar, StellarMasses.metal, StellarMasses.absMag_r<br>\n",
      "FROM SpecObj JOIN StellarMasses ON SpecObj.CATAID = StellarMasses.CATAID<br>\n",
      "WHERE  SpecObj.nq>2 and SpecObj.prob >0.75<br>\n",
      "\n",
      "The purpose of the JOIN syntax in the FROM clause is to make sure we are looking at the same galaxies in each table. The conditions in the WHERE clause make sure that the galaxies spectrum was high quality and that the redshift was certain.\n",
      "\n",
      "Using the <a href=\"http://astroquery.readthedocs.org/en/latest/#\">astroquery module</a> we can query the GAMA database inside python and return the results in a data table (which in the astronomical tradition of clever nameing, I call dataTable)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from astroquery.gama import GAMA\n",
      "from astropy.table import Table,Column\n",
      "from astropy import units as u\n",
      "from astropy.coordinates import ICRS, CartesianRepresentation, SkyCoord, Distance\n",
      "#Acquire Data from GAMA database\n",
      "dataTable = GAMA.query_sql('\\\n",
      "   SELECT  SpecObj.ra, SpecObj.dec, SpecObj.z,\\\n",
      "          StellarMasses.logmstar, StellarMasses.metal, StellarMasses.absMag_r\\\n",
      "   FROM SpecObj JOIN StellarMasses ON SpecObj.CATAID = StellarMasses.CATAID\\\n",
      "   WHERE  SpecObj.nq>2 and SpecObj.prob >0.75'\\\n",
      ")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Downloading http://www.gama-survey.org/dr2/query/../tmp/GAMA_5TARzH.fits"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " [Done]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Step #2 Calculate distance to the galaxies from their redshift"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To plot the galaxies in 3D space we need both the position on the sky and a distance. We will convert the redshift into a distance, but how? On cosmological scales in an expanding universe with a finite speed of light, there are several ways to measure distance. Which one should we use? People have made different choices here. Some use the lookback time distance, which is the light travel time from the galaxy multiplied by the speed of light. Others have used the luminosity distance which is the distance at which the galaxy would have its observed brightness given the inverse square law. I believe that the best choice is the co-moving distance. This distance is the one that best matches our expectation of what distance is. If we were able to stop the expansion of the universe and go and measure the distance to galaxies that is what we would get.\n",
      "\n",
      "Yes, that does mean that you will get distances greater than 13.8 billion light years, but that is OK, really it is. If you don't believe me see the Wikipedia page about misconceptions regarding the size of the Universe.\n",
      "\n",
      "Now to calculate the comoving distance we'll need to choose a set of Cosmological parameters. One of the really nice features of astropy is  the built in set of standard cosmologies. Here I choose the Planck13 cosmology defined <a href=\"http://arxiv.org/abs/1303.5076\">here</a>. It just takes three lines of code to choose the cosmological parameters, calculate the distances and add a new column to the data table.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from astropy.cosmology import Planck13\n",
      "#Calculate a new column of comoving distance to the galaxies\n",
      "distCol=Column(Planck13.comoving_distance(dataTable['z']),name='comoving_distance',meta={'ucd': 'pos.distance'})\n",
      "#Add that column to the data table (as the third column after ra and dec)\n",
      "dataTable.add_column(distCol,2)\n",
      "#Create a column of coordinate objects because it will be useful later\n",
      "coordsCol=SkyCoord(dataTable['ra'],dataTable['DEC'],unit=(u.degree, u.degree),\\\n",
      "                   distance=Distance(dataTable['comoving_distance'],u.Mpc),frame='icrs')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Print the Catalog"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print len(dataTable)\n",
      "dataTable"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "65887\n"
       ]
      },
      {
       "html": [
        "<table id=\"table61473008\"><thead><tr><th>ra</th><th>DEC</th><th>comoving_distance</th><th>z</th><th>logmstar</th><th>metal</th><th>absMag_r</th></tr></thead><thead><tr><th></th><th></th><th>Mpc</th><th></th><th></th><th></th><th></th></tr></thead><tr><td>174.006</td><td>0.72093</td><td>220.931793244</td><td>0.05054</td><td>9.12851042187</td><td>0.0330536359134</td><td>-18.1912828065</td></tr><tr><td>174.02279</td><td>0.70594</td><td>1348.23828782</td><td>0.33119</td><td>11.0503601715</td><td>0.0202665491508</td><td>-22.8990099392</td></tr><tr><td>174.10071</td><td>0.65891</td><td>960.638445368</td><td>0.22979</td><td>10.847068929</td><td>0.0184374417908</td><td>-21.7511866961</td></tr><tr><td>174.15312</td><td>0.81543</td><td>16.574370445</td><td>0.00375</td><td>7.48716971791</td><td>0.000835064700905</td><td>-15.5538521806</td></tr><tr><td>174.18493</td><td>0.70906</td><td>316.458111207</td><td>0.07278</td><td>10.1603001507</td><td>0.0153691129223</td><td>-20.5006711693</td></tr><tr><td>174.2805</td><td>0.70608</td><td>490.601220644</td><td>0.11397</td><td>9.55955138127</td><td>0.0130184432258</td><td>-19.7482673646</td></tr><tr><td>174.30275</td><td>0.78998</td><td>322.692558895</td><td>0.07424</td><td>9.82172906009</td><td>0.0066186863944</td><td>-20.0559795935</td></tr><tr><td>174.30554</td><td>0.79034</td><td>323.930364114</td><td>0.07453</td><td>9.47458337058</td><td>0.00753276230635</td><td>-19.4336258665</td></tr><tr><td>174.33537</td><td>0.83656</td><td>460.129001741</td><td>0.1067</td><td>9.85374077052</td><td>0.0148749711966</td><td>-20.1791723545</td></tr><tr><td>174.34688</td><td>0.69645</td><td>814.85045396</td><td>0.19309</td><td>10.7435902879</td><td>0.0116710343688</td><td>-21.785559543</td></tr><tr><td>174.39605</td><td>0.82075</td><td>420.898291295</td><td>0.09738</td><td>10.5514539145</td><td>0.0109531939707</td><td>-21.1061687293</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>221.0395</td><td>0.7408</td><td>1089.99846388</td><td>0.26299</td><td>10.6328220566</td><td>0.0132416501559</td><td>-21.6004346227</td></tr><tr><td>221.35471</td><td>0.66662</td><td>0.309652244074</td><td>7e-05</td><td>-3.78759898628</td><td>0.00513032257297</td><td>13.0997089078</td></tr><tr><td>221.41271</td><td>0.79681</td><td>0.840461319196</td><td>0.00019</td><td>-3.22269328729</td><td>0.0123144631116</td><td>13.5748862713</td></tr><tr><td>221.41892</td><td>0.80136</td><td>694.724548139</td><td>0.16339</td><td>10.4246441576</td><td>0.0171651008519</td><td>-20.625278007</td></tr><tr><td>222.70867</td><td>0.77829</td><td>2.25580783759</td><td>0.00051</td><td>-1.42445525549</td><td>0.00589734235856</td><td>7.91604008527</td></tr><tr><td>182.92743</td><td>1.06032</td><td>1206.99497028</td><td>0.29356</td><td>9.98784676898</td><td>0.0203818348799</td><td>-22.7466977624</td></tr><tr><td>183.46533</td><td>1.21717</td><td>1.45970129094</td><td>0.00033</td><td>-2.89053379897</td><td>0.0135255803682</td><td>13.1098955275</td></tr><tr><td>211.55021</td><td>1.10102</td><td>4.20156704942</td><td>0.00095</td><td>5.38068476844</td><td>0.0146118260433</td><td>-6.66670411714</td></tr><tr><td>211.61046</td><td>1.05388</td><td>4.99744494288</td><td>0.00113</td><td>5.76503529501</td><td>0.00617013977633</td><td>-9.3358116772</td></tr><tr><td>212.0575</td><td>1.0527</td><td>327.856048501</td><td>0.07545</td><td>9.57088001257</td><td>0.00836464776078</td><td>-19.0808277015</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "<Table rows=65887 names=('ra','DEC','comoving_distance','z','logmstar','metal','absMag_r') units=(None,None,'Mpc',None,None,None,None)>\n",
        "array([ (174.006, 0.7209299802780151, 220.93179324389214, 0.050540000200271606, 9.12851042187451, 0.0330536359133739, -18.1912828065272),\n",
        "       (174.02279, 0.7059400081634521, 1348.238287816697, 0.33118999004364014, 11.0503601714543, 0.020266549150793, -22.8990099391754),\n",
        "       (174.10071, 0.6589099764823914, 960.6384453680855, 0.22979000210762024, 10.8470689289949, 0.0184374417907734, -21.7511866961371),\n",
        "       ...,\n",
        "       (211.55021, 1.1010199785232544, 4.201567049423963, 0.0009500000160187483, 5.38068476844009, 0.0146118260433305, -6.66670411713684),\n",
        "       (211.61046, 1.053879976272583, 4.99744494287576, 0.001129999989643693, 5.76503529501127, 0.00617013977633484, -9.33581167720381),\n",
        "       (212.0575, 1.0527000427246094, 327.85604850058616, 0.07545000314712524, 9.57088001256696, 0.0083646477607815, -19.0808277014637)], \n",
        "      dtype=[('ra', '>f8'), ('DEC', '>f4'), ('comoving_distance', '<f8'), ('z', '>f4'), ('logmstar', '>f8'), ('metal', '>f8'), ('absMag_r', '>f8')])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that we have calculated the distances we are ready to visualize the dataset. The process will vary a bit depending on which software package we are using. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Visualizing the dataset in Microsoft Research's Worldwide Telescope"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can then send the data table directly to WWT as a new layer utilizing WWT's layer control API (<a href=\"http://www.worldwidetelescope.org/Docs/WorldWideTelescope_lcapi.html\">LCAPI</a>). It can be sent to a WWT client running on the same machine or on a remote machine. For example I can run this notebook on my laptop and send data to Adler's planetarium dome. John Zuhone's <a href=\"http://www.jzuhone.com/pywwt/index.html\">pyWWT</a> package provides a python interface to the LCAPI as well as some nice utilities (for example one to help create colormaps)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pywwt.mods import * #import the pyWWT package\n",
      "#Connect to WWT\n",
      "my_wwt = WWTClient(host=\"127.0.0.1\") #Make sure that remote connections are enabled on the WWT client"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Once the connection is made all we have to do is set up the new layer, set the visualization properties for the layer, and import out data table into the layer."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Set up WWT layer\n",
      "GAMA_layer = my_wwt.new_layer(\"Sky\", \"GAMA Galaxies\", dataTable.colnames)\n",
      "#Set visualization parameters in WWT\n",
      "props_dict = {\"CoordinatesType\":\"Spherical\",\\\n",
      "              \"MarkerScale\":\"Screen\",\\\n",
      "              \"PointScaleType\":\"Constant\",\\\n",
      "              \"ScaleFactor\":\"2\",\\\n",
      "              \"ShowFarSide\":\"True\",\\\n",
      "              \"AltColumn\":\"2\",\\\n",
      "              \"AltType\":\"Distance\",\\\n",
      "              \"AltUnit\":\"MegaParsecs\",\\\n",
      "              \"RaUnits\":\"Degrees\",\\\n",
      "              \"TimeSeries\":\"False\"}\n",
      "GAMA_layer.set_properties(props_dict)\n",
      "#Send data to WWT client\n",
      "GAMA_layer.update(data=dataTable, purge_all=True, no_purge=False, show=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The result looks like this:\n",
      "   <img src=\"https://raw.githubusercontent.com/IPSScienceVisualization/python-tutorials/master/images/WWT2.PNG\"><br>\n",
      "Switch to Sky Mode, change the marker type to circles and scale them up a bit and and you have this:<img src=\"https://raw.githubusercontent.com/IPSScienceVisualization/python-tutorials/master/images/WWT3.PNG\"><br>"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}