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",
    "# For downloading files\n",
    "from astropy.utils.data import download_file\n",
    "from astropy.io import fits\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 result as an Astropy Table makes it easier to look at the contents, to call the service itself, we cannot use the row of that table.  You have to use the entry in the service result list itself.  So use the table to browse, but select the list of services itself using the properties that have been defined as attributes such as short_name and ivoid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galex_stsci=[s for s in uv_services if 'GALEX' in s.short_name and 'stsci' in s.ivoid][0]\n",
    "galex_heasarc=[s for s in uv_services if 'GALEX' in s.short_name and 'heasarc' in s.ivoid][0]"
   ]
  },
  {
   "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.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": [
    "## If you only run this once, you can do it in memory in one line:\n",
    "##  This fetches the FITS as an astropy.io.fits object in memory\n",
    "#dataobj=im_table_heasarc[0].getdataobj()\n",
    "## But if you might run this notebook repeatedly with limited bandwidth, \n",
    "##  download it once and cache it.  \n",
    "file_name = download_file(im_table_heasarc[0].getdataurl(),cache=True)\n",
    "dataobj=fits.open(file_name)\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",
    "## See above for options\n",
    "#xhdu_list=first_fits_image_row.getdataobj()\n",
    "file_name = download_file(first_fits_image_row.getdataurl(),cache=True)\n",
    "xhdu_list=fits.open(file_name)\n",
    "\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 and look at it:\n",
    "#hdu_list=spec_tables[0].getdataobj()\n",
    "file_name = download_file(spec_tables[0].getdataurl(),cache=True)\n",
    "hdu_list=fits.open(file_name)\n",
    "\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
}