Science Use Case - Proposal Prep pyvo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing a proposal\n",
    "\n",
    "**<font color=red>Updated to use [pyvo pulled from GitHub](https://github.com/astropy/pyvo) on 2019-07-08.  Notes in red about issues or features we can/should add to pyvo from our navo_utils.</font>**\n",
    "\n",
    "The Story: Suppose that you are preparing to write a proposal on NGC1365, aiming to investigate the intriguing black hole spin this galaxy with Chandra grating observations (see: https://www.space.com/19980-monster-black-hole-spin-discovery.html ) \n",
    "\n",
    "In writing proposals, there are often the same tasks that are required: including finding and analyzing previous observations of the proposal, and creating figures that include, e.g., multiwavelength images and spectrum for the source. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "##  For DEBUGGING\n",
    "## \n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# As a hint, we include the code block for Python modules that you will likely need to import:   \n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline  \n",
    "\n",
    "##  Use pyvo as much as possible\n",
    "import pyvo as vo\n",
    "\n",
    "\n",
    "## There are a number of relatively unimportant warnings that \n",
    "## show up, so for now, suppress them:\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Find out what the previously quoted Chandra 2-10 keV flux of the central source is for NGC 1365.  \n",
    "\n",
    "Hint: Do a table access protocol (TAP) search in the heasarc, to find potential table with this information\n",
    "\n",
    "<font color=red>No \"source\" but \"keywords\" is fine for the same thing.  \n",
    "    \n",
    "Then you can take the set of results and browse them as a table with \n",
    "\n",
    "    tap_services.table\n",
    "    \n",
    "or you can take one of them and use it as a service with\n",
    "\n",
    "    tap_services[0].search(query)\n",
    "\n",
    "(Is there another way to get the URL as a string?  The getdataurl() doesn't work with this object.  I guess the idea is that you shouldn't *need* to see the URL.)  \n",
    "</font>\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#tap_services=Registry.query(service_type='table',source='heasarc')\n",
    "tap_services=vo.regsearch(servicetype='table',keywords=['heasarc'])\n",
    "url=tap_services[0]['access_url'].decode(\"utf-8\")\n",
    "print(url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: The Chansngcat ( https://heasarc.gsfc.nasa.gov/W3Browse/chandra/chansngcat.html ) table is likely the best table.  Create a table with ra, dec, exposure time, and flux (and flux errors) from the public.chansngcat catalog for Chandra observations matched within 0.1 degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Found {} services:\".format(len(tap_services)))\n",
    "for s in tap_services:\n",
    "    print(s.res_description)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red>\n",
    "I think we should be learning/teaching how to chain pyvo objects themselves.  For example, ** \n",
    "\n",
    "    tap_services[0]\n",
    "    \n",
    "is a row of the registry results, and if you tack on .service, it becomes the appropriate service type object (e.g., a TAPService or an SIAService).  But you don't necessarily have to, because special functions work on the registry result row itself.  These two are equivalent I think:\n",
    "\n",
    "    tap_services[0].service.run_sync(query)\n",
    "    tap_services[0].search(query)\n",
    "    \n",
    "for example for a TAP service, or \n",
    "\n",
    "    image_services[0].search(position,size)\n",
    "\n",
    "for an image.\n",
    "\n",
    "One annoying exception that I think should be \"exposed\" higher is that for a TAP services, if you want to find the available tables, then that's the .tables object, but this doesn't work on a row of a result set;  you have to specify the .service attribute first:\n",
    "\n",
    "    tap_services[0].service.tables\n",
    "    \n",
    "(not to be confused with .table, which converts the results into an Astropy Table object) which returns an VOSITables object that can work like a dictionary whose keys are the table names and values are tables of its meta data.  So you can do the following, or at least should be able to when it's made to work(?):\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heasarc_tables=tap_services[0].service.tables\n",
    "for tablename in heasarc_tables.keys():\n",
    "    if \"chan\" in tablename:  \n",
    "        print(tablename)\n",
    "        print(heasarc_tables[tablename].description)\n",
    "        ## These don't appear to work\n",
    "        #heasarc_tables[tablename].describe()\n",
    "        #print(heasarc_tables[tablename].columns)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red>  But is this starting to be a bit annoying for the average user, making sure to use the right methods to get the right kinds of objects called in the right ways?  I've figured this out reading the code and lots of trial and error.  The docs don't cover this quite the way I needed.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the coordinate for NGC 1365\n",
    "import astropy.coordinates as coord\n",
    "pos=coord.SkyCoord.from_name(\"ngc1365\")\n",
    "query=\"\"\"SELECT ra, dec, exposure, flux, flux_lower, flux_upper \n",
    "    FROM public.chansngcat as cat \n",
    "    where \n",
    "    contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{},{},0.1))=1 \n",
    "    and \n",
    "    cat.exposure > 0 order by cat.exposure\n",
    "    \"\"\".format(pos.ra.deg, pos.dec.deg)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Show the results from your table query for NGC 1365 \n",
    "#results=tap_services[0].service.run_sync(query)\n",
    "results=tap_services[0].search(query)\n",
    "results.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Step 2: Make Images: \n",
    "\n",
    "#### Create ultraviolet and X-ray images\n",
    "Hint: Start by checking what UV image services exist (e.g., GALEX?)\n",
    "\n",
    "<font color=red>\n",
    "As above, add .table to browse the results more easily.\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Note that to browse more 'prettily', use the .table attribute\n",
    "uv_services=vo.regsearch(servicetype='image',keywords='galex', waveband='uv')\n",
    "## Can't do two columns at a time with this way:\n",
    "uv_services.table['ivoid','short_name']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how you select the services you're interested in.  Take the table of results\n",
    "\n",
    "    uv_services.table\n",
    "\n",
    "and search the relevant column for the relevant string\n",
    "\n",
    "    np.isin(uv_services.table['short_name'],b'GALEX')\n",
    "\n",
    "and use the resulting indices of where the GALEX entries are to look at three of the table columns for only those GALEX rows:\n",
    "\n",
    "    uv_services.table[\n",
    "        np.isin(uv_services.table['short_name'],b'GALEX')\n",
    "        ]['ivoid','short_name','access_url']\n",
    "\n",
    "\n",
    "There are two that serve GALEX images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uv_services.table[\n",
    "    np.isin(uv_services.table['short_name'],b'GALEX')\n",
    "    ]['ivoid','short_name','access_url']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red>But the above uses and returns a table which you cannot call.  So you want to get the elements of the original object that are callable.\n",
    "\n",
    "But slicing the object as though it's a table doesn't work, i.e., these don't work:\n",
    "\n",
    "    services[0:1]\n",
    "    services[ np.where(np.isin(services.table['short_name'],b'GALEX'))[0] ]\n",
    "\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## The where returns a tuple, but we only want the first part of the tuple, the array of indices.\n",
    "##  But you cannot pass more than one integer to a results object to get something callable, i.e.\n",
    "##  though some table operations work, it doesn't slice like a table.  So get them separately by\n",
    "##  specifying the first and second of the two matches from the where:\n",
    "galex_stsci=uv_services[int(np.where(np.isin(uv_services.table['short_name'],b'GALEX'))[0][0])]\n",
    "galex_heasarc=uv_services[int(np.where(np.isin(uv_services.table['short_name'],b'GALEX'))[0][1])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font color=red>But this is overly ugly.  Took me lots of trial and error and reading code to figure out how exactly to get what I want out of the pyvo objects.  And given that the results may be returned in a different order tomorrow, we cannot just hard-wire the numbers.  There's got to be a better way to do this.</font>**\n",
    "\n",
    "<font color=red>  This is how we did it with the navo_utils:\n",
    "\n",
    "    uv_services=Registry.query(service_type='image',keyword='galex', waveband='uv')\n",
    "    im_gal = uv_services[np.isin(uv_services['short_name'], \n",
    "                                        ['GALEX'])]\n",
    "    \n",
    "Since everything was just tables, they could be sliced and diced.  But then you cannot do the below way of just tacking on .search() to the result entry, which is handy.\n",
    "</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: Next create a UV image for the source "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do an image search for NGC 1365 in the UV service found above\n",
    "im_table_stsci=galex_stsci.search(pos=pos,size=0.1)\n",
    "im_table_stsci.table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do an image search for NGC 1365 in the UV service found above\n",
    "im_table_heasarc=galex_heasarc.search(pos=pos,size=0.2)\n",
    "im_table_heasarc.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font color=red>Note that our navo_utils.Image.query() takes a position and RADIUS.  This takes a position and a SIZE, which is apparently a diameter.  So there's a factor of two difference in what you give.  Is this what people will expect from an SIA interface?</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Get the FITS file, first of above rows.  \n",
    "\n",
    "##  There are a couple of things you can extract \n",
    "##   from each row of the results:\n",
    "\n",
    "\n",
    "## Shouldn't need this\n",
    "print(im_table_heasarc[0].getdataurl())  \n",
    "## For downloading the file, see below\n",
    "dataset=im_table_heasarc[0].getdataset() \n",
    "## For grabbing the FITS as an astropy.io.fits object in memory\n",
    "dataobj=im_table_heasarc[0].getdataobj()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(dataset)) \n",
    "print(type(dataobj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make an image: \n",
    "from pylab import figure, cm\n",
    "from matplotlib.colors import LogNorm\n",
    "plt.matshow(dataobj[0].data, origin='lower', cmap=cm.gray_r, norm=LogNorm(vmin=0.005, vmax=0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font color=red>Note that each record has a lot of information in it.  Need a method to print this information in a useful way, which I don't currently see.  Doesn't have, e.g., a .table attribute.</font>**\n",
    "\n",
    "**<font color=red>The change in the column names (e.g., \"Format\" versus \"content-Type\" versus \"imgfmt\" etc.) is presumably why there's a getdataformat() function.  Not sure I like this, though.  Can see how it relates to our navo_utils UCD conversion stuff.</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(im_table_heasarc[0]))\n",
    "print(im_table_heasarc[0].getdataformat())\n",
    "dir(im_table_heasarc[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: Repeat steps for X-ray image. (Note: Ideally, we would find an image in the Chandra 'cxc' catalog) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_services=vo.regsearch(servicetype='image',keywords=['chandra'], waveband='x-ray')\n",
    "print(x_services.table['short_name','ivoid'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Do an image search for NGC 1365 in the X-ray CDA service found above\n",
    "xim_table=x_services[0].search(pos=pos,size=0.2)\n",
    "## Some of these are FITS and some JPEG.  Look at the columns:\n",
    "print( xim_table.table.columns )\n",
    "first_fits_image_row=xim_table[ int(np.where(np.isin(xim_table.table['imgfmt'],b'image/fits'))[0][0]) ] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an image from the first FITS file (index=1) by downloading:\n",
    "xhdu_list=first_fits_image_row.getdataobj()\n",
    "\n",
    "plt.imshow(xhdu_list[0].data, origin='lower', cmap='cool', norm=LogNorm(vmin=0.1, vmax=500.))\n",
    "plt.xlim(460, 560)\n",
    "plt.ylim(460, 560)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Step 3: Make a spectrum: \n",
    "\n",
    "#### Find what Chandra spectral observations exist already for this source. \n",
    "Hint: try searching for X-ray spectral data tables using the registry query"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font color=red>Why do these two ways of asking for spectral access services return different results?  Looks like it's doing an exact rather than partial match, so services that are \"x-ray#gamma-ray\" don't get found by pyvo.  Didn't help to try waveband='%x-ray%', but there ought to be a way.</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the TAP protocol to list services that contain X-ray spectral data\n",
    "xsp_services=vo.regsearch(servicetype='ssa',waveband='x-ray')\n",
    "xsp_services.table['short_name','ivoid','waveband']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red>This is the way our navo_utils.registry did it and what it finds:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from navo_utils.registry import Registry\n",
    "services=Registry.query(service_type='spectr',waveband='x-ray')\n",
    "services['short_name','ivoid','waveband']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint 2: Take a look at what data exist for our candidate, NGC 1365."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font color=red>I'd also request that \"service unavailable\" and a few other errors be trapped for and a simple message printed rather than the whole traceback from an exception</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_tables=xsp_services[0].search(pos=pos,radius=0.2,verbose=True)\n",
    "spec_tables.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint 3: Download the data to make a spectrum. Note: you might end here and use Xspec to plot and model the spectrum. Or ... you can also try to take a quick look at the spectrum. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font color=red>This ends up with a file named \"acisf13920N001_pha2.None\" instead of the correct \".fits\" suffix.  And ought to be wrapped as something simpler so the user doesn't have to know about .make_dataset_filename() and .getdataset().read() etc.</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Or write it to disk\n",
    "with open('downloads/{}'.format(spec_tables[0].make_dataset_filename()),'wb') as outfile:\n",
    "    outfile.write(spec_tables[0].getdataset().read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu_list=spec_tables[0].getdataobj()\n",
    "spectra=hdu_list[1].data\n",
    "print(spectra.columns)\n",
    "print(len(spectra))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extension: Making a \"quick look\" spectrum. For our purposes, the 1st order of the HEG grating data would be sufficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "j=1\n",
    "for i in range(len(spectra)):\n",
    "    matplotlib.rcParams['figure.figsize'] = (8, 3)\n",
    "    if abs(spectra['TG_M'][i]) == 1 and (spectra['TG_PART'][i]) == 1:\n",
    "        ax=plt.subplot(1,2,j)\n",
    "        pha = plt.plot( spectra['CHANNEL'][i],spectra['COUNTS'][i])\n",
    "        ax.set_yscale('log')\n",
    "        if spectra['TG_PART'][i] == 1:\n",
    "            instr='HEG'\n",
    "        ax.set_title(\"{grating}{order:+d}\".format(grating=instr, order=spectra['TG_M'][i]))\n",
    "        plt.tight_layout()\n",
    "        j=j+1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html).  (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.ipynb) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations! You have completed this notebook exercise."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.7"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}