SOLUTIONS-ProposalPrep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing a proposal\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": [
    "# 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",
    "import pyvo as vo\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\", module=\"astropy.io.votable.*\")\n",
    "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  This gets all services matching, which should be a list of only one:\n",
    "tap_services=vo.regsearch(servicetype='table',keywords=['heasarc'])\n",
    "#  This fetches the list of tables that this service serves:\n",
    "heasarc_tables=tap_services[0].service.tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for tablename in heasarc_tables.keys():\n",
    "    if \"chan\" in tablename:  \n",
    "        print(\"Table {} has columns={}\\n\".format(\n",
    "            tablename,\n",
    "            sorted([k.name for k in heasarc_tables[tablename].columns ])))"
   ]
  },
  {
   "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": [
    "# Get the coordinate for NGC 1365\n",
    "import astropy.coordinates as coord\n",
    "pos=coord.SkyCoord.from_name(\"ngc1365\")\n",
    "# Construct a query that will get the ra, dec, exposure time, flux, and flux errors \n",
    "#  from this catalog in the region around this source:\n",
    "query=\"\"\"SELECT ra, dec, exposure, flux, flux_lower, flux_upper FROM public.chansngcat as cat \n",
    "    where contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{},{},0.1))=1 \n",
    "    and cat.exposure > 0 order by cat.exposure\"\"\".format(pos.ra.deg, pos.dec.deg)\n",
    "#  Submit the query.  (See the CS_Catalog_queries.ipynb for\n",
    "#    information about these two search options.)\n",
    "#results=tap_services[0].service.run_async(query)\n",
    "results=tap_services[0].search(query)\n",
    "#  Look at the results\n",
    "results.to_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?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Note that to browse the columns, use the .table attribute\n",
    "uv_services=vo.regsearch(servicetype='image',keywords='galex', waveband='uv')\n",
    "uv_services.to_table()['ivoid','short_name']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The keyword search for 'galex' returned a bunch of things that may have mentioned it, but let's just use the ones that have GALEX as their short name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uv_services.to_table()[\n",
    "    np.isin(uv_services.to_table()['short_name'],b'GALEX')\n",
    "    ]['ivoid','short_name','access_url']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though using the .table attribute makes it easier to look at the contents, to call the service itself, we cannot use the row of that table but the row of the service result list itself.  So use the table to find the one you want, then use that as the index of the list of services:"
   ]
  },
  {
   "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.to_table()['short_name'],b'GALEX'))[0][0])]\n",
    "galex_heasarc=uv_services[int(np.where(np.isin(uv_services.to_table()['short_name'],b'GALEX'))[0][1])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: Next create a UV image for the source "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "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.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#  Let's see what HEASARC offers, and this time limit it to FITS \n",
    "#   this option doesn't currently work for STScI's service)\n",
    "im_table_heasarc=galex_heasarc.search(pos=pos,size=0.1,format='image/fits')\n",
    "im_table_heasarc.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## For grabbing the FITS as an astropy.io.fits object in memory\n",
    "dataobj=im_table_heasarc[0].getdataobj()\n",
    "print(type(dataobj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the FITS file (which is index 0 for the NUV image or index=2 for the FUV 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": [
    "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.to_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.to_table().columns )\n",
    "first_fits_image_row=xim_table[ int(np.where(np.isin(xim_table.to_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": "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.to_table()['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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_tables=xsp_services[0].search(pos=pos,radius=0.2,verbose=True)\n",
    "spec_tables.to_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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Get it in memory and look at it:\n",
    "hdu_list=spec_tables[0].getdataobj()\n",
    "spectra=hdu_list[1].data\n",
    "print(spectra.columns)\n",
    "print(len(spectra))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Or write it to disk\n",
    "import os\n",
    "if not os.path.isdir('downloads'):\n",
    "    os.makedirs(\"downloads\")\n",
    "fname=spec_tables[0].make_dataset_filename()\n",
    "#  Known issue where the suffix is incorrect:\n",
    "fname=fname.replace('None','fits')\n",
    "with open('downloads/{}'.format(fname),'wb') as outfile:\n",
    "    outfile.write(spec_tables[0].getdataset().read())"
   ]
  },
  {
   "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": 4
}