SOLUTIONS-InspectingACandidateList.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Science User Case - Inspecting a Candidate List\n",
    "\n",
    "Ogle et al. (2016) mined the NASA/IPAC Extragalactic Database (NED) to identify a new type of galaxy: Superluminous Spiral Galaxies.\n",
    "\n",
    "Here's the paper: https://ui.adsabs.harvard.edu//#abs/2016ApJ...817..109O/abstract\n",
    "\n",
    "Table 1 lists the positions of these Super Spirals. Based on those positions, let's create multiwavelength cutouts for each super spiral to see what is unique about this new class of objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Import the Python modules we'll be using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Suppress unimportant warnings.\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n",
    "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")\n",
    "warnings.filterwarnings('ignore', '.*RADECSYS=*', append=True)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# For downloading files\n",
    "from astropy.utils.data import download_file\n",
    "\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.io import fits\n",
    "from astropy.nddata import Cutout2D\n",
    "import astropy.visualization as vis\n",
    "from astropy.wcs import WCS\n",
    "from astroquery.ned import Ned\n",
    "\n",
    "import pyvo as vo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell prepares the notebook to display our visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Search NED for objects in this paper.\n",
    "\n",
    "Insert a Code Cell below by clicking on the \"Insert\" Menu and choosing \"Insert Cell Below\". Then consult QuickReference.md to figure out how to use astroquery to search NED for all objects in a paper, based on the refcode of the paper. Inspect the resulting astropy table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objects_in_paper = Ned.query_refcode('2016ApJ...817..109O')\n",
    "objects_in_paper.show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Filter the NED results.\n",
    "\n",
    "The results from NED will include galaxies, but also other kinds of objects. Print the 'Type' column to see the full range of classifications. Next, print the 'Type' of just the first source in the table, in order to determine its data type (since Python 3 distinguishes between strings and byte strings). Finally, use the data type information to filter the results so that we only keep the galaxies in the list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objects_in_paper['Type']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's see whether we are dealing with strings or byte strings\n",
    "objects_in_paper[0]['Type']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Byte strings it is!\n",
    "galaxies = objects_in_paper[objects_in_paper['Type'] == b'G']\n",
    "\n",
    "galaxies.show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Search the NAVO Registry for image resources.\n",
    "\n",
    "The paper selected super spirals using WISE, SDSS, and GALEX images. Search the NAVO registry for all image resources, using the 'service_type' search parameter. How many image resources are currently available?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_services = vo.regsearch(servicetype='image')\n",
    "\n",
    "print(f'{len(image_services)} result(s) found.')\n",
    "\n",
    "image_services.to_table()['ivoid', 'short_name', 'res_title']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Search the NAVO Registry for image resources that will allow you to search for AllWISE images.\n",
    "\n",
    "There are hundreds of image resources...too many to quickly read through. Try adding the 'keywords' search parameter to your registry search, and find the image resource you would need to search the AllWISE images. Remember from the Known Issues that 'keywords' must be a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allwise_image_services = vo.regsearch(servicetype='image', keywords=['allwise'])\n",
    "\n",
    "print(f'{len(allwise_image_services)} result(s) found.')\n",
    "\n",
    "allwise_image_services.to_table()['ivoid', 'short_name', 'res_title']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Choose the AllWISE image service that you are interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allwise_image_service = allwise_image_services[0]\n",
    "allwise_image_service.service"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Choose one of the galaxies in the NED list.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra = galaxies['RA'][0]\n",
    "dec = galaxies['DEC'][0]\n",
    "pos = SkyCoord(ra, dec, unit = 'deg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Search for a list of AllWISE images that cover this galaxy.\n",
    "\n",
    "How many images are returned? Which are you most interested in?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allwise_image_table = allwise_image_service.search(pos=pos, size=0)\n",
    "allwise_image_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Use the .to_table() method to view the results as an Astropy table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allwise_images = allwise_image_table.to_table()\n",
    "allwise_images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. From the result in 8., select the first record for an image taken in WISE band W1 (3.6 micron)\n",
    "\n",
    "Hints:\n",
    "* Loop over records and test on the `.bandpass_id` attribute of each record\n",
    "* Print the `.title` and `.bandpass_id` of the record you find, to verify it is the right one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for allwise_image_record in allwise_image_table:\n",
    "    if 'W1' in allwise_image_record.bandpass_id:\n",
    "        break\n",
    "print(allwise_image_record.title, allwise_image_record.bandpass_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Visualize this AllWISE image."
   ]
  },
  {
   "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",
    "#allwise_w1_image = allwise_image_record.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(allwise_image_record.getdataurl(), cache=True)  \n",
    "allwise_w1_image = fits.open(file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "wcs = WCS(allwise_w1_image[0].header)\n",
    "ax = fig.add_subplot(1, 1, 1, projection=wcs)\n",
    "ax.imshow(allwise_w1_image[0].data, cmap='gray_r', origin='lower', vmax = 10)\n",
    "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Plot a cutout of the AllWISE image, centered on your position.\n",
    "\n",
    "Try a 60 arcsecond cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 60\n",
    "cutout = Cutout2D(allwise_w1_image[0].data, pos, (size, size), wcs=wcs)\n",
    "wcs = cutout.wcs\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(1, 1, 1, projection=wcs)\n",
    "ax.imshow(cutout.data, cmap='gray_r', origin='lower', vmax = 10)\n",
    "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13. Try visualizing a cutout of a GALEX image that covers your position.\n",
    "\n",
    "Repeat steps 4, 5, 6, 8 through 12 for GALEX."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galex_image_services = vo.regsearch(keywords=['galex'], servicetype='image')\n",
    "print(f'{len(galex_image_services)} result(s) found.')\n",
    "galex_image_services.to_table()['ivoid', 'short_name', 'res_title']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galex_image_service = galex_image_services[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galex_image_table = galex_image_service.search(pos=pos, size=0.0, intersect='covers')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(galex_image_table)):\n",
    "    if (('image/fits' in galex_image_table[i].format) and\n",
    "        (galex_image_table['enrValue'][i]==2.35e-07)):\n",
    "        break\n",
    "galex_image_record = galex_image_table[i]\n",
    "print(galex_image_record.title, galex_image_record.bandpass_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## See above regarding two ways to do this: \n",
    "#galex_nuv_image = fits.open(galex_image_record.getdataurl())\n",
    "file_name = download_file(galex_image_record.getdataurl(), cache=True)  \n",
    "galex_nuv_image=fits.open(file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_data = galex_nuv_image[0].data\n",
    "print('Min:', np.min(image_data))\n",
    "print('Max:', np.max(image_data))\n",
    "print('Mean:', np.mean(image_data))\n",
    "print('Stdev:', np.std(image_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1, projection=WCS(galex_nuv_image[0].header))\n",
    "ax.imshow(galex_nuv_image[0].data, cmap='gray_r', origin='lower', vmin=0.0, vmax=0.01)\n",
    "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutout = Cutout2D(galex_nuv_image[0].data, pos, size, wcs=WCS(galex_nuv_image[0].header))\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(1, 1, 1, projection=cutout.wcs)\n",
    "ax.imshow(cutout.data, cmap='gray_r', origin='lower', vmin = 0.0, vmax = 0.01)\n",
    "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 14. Try visualizing a cutout of an SDSS image that covers your position.\n",
    "\n",
    "Hints:\n",
    "* Search the registry using `keywords=['sloan']\n",
    "* Find the service with a `short_name` of `b'SDSS SIAP'`\n",
    "* From Known Issues, recall that an empty string must be specified to the `format` parameter dues to a bug in the service.\n",
    "* After obtaining your search results, select r-band images using the `.title` attribute of the records that are returned, since `.bandpass_id` is not populated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_image_services = vo.regsearch(keywords=['sloan'], servicetype='image')\n",
    "sdss_image_services.to_table()['ivoid', 'short_name', 'res_title', 'source_value']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Use list comprehension to check each service's short_name attribute.\n",
    "#   Given the above, we know the first match is the right one.  \n",
    "sdss_image_service = [s for s in sdss_image_services if 'SDSS SIAP' in s.short_name ][0]\n",
    "sdss_image_service.short_name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a workaround to a bug in the SDSS service, pass `format=''` when searching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_image_table = sdss_image_service.search(pos=pos, size=0.0, format='', intersect='covers')\n",
    "len(sdss_image_table['Title'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sdss_rband_record in sdss_image_table:\n",
    "    if 'Sloan Digital Sky Survey - Filter r' in sdss_rband_record.title:\n",
    "        break\n",
    "print(sdss_rband_record.title, sdss_rband_record.bandpass_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##  See above regarding two ways to do this\n",
    "# sdss_rband_image = fits.open(sdss_rband_record.getdataurl())\n",
    "file_name = download_file(sdss_rband_record.getdataurl(), cache=True)  \n",
    "sdss_rband_image=fits.open(file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(1, 1, 1, projection=WCS(sdss_rband_image[0].header))\n",
    "\n",
    "interval = vis.PercentileInterval(99.9)\n",
    "vmin,vmax = interval.get_limits(sdss_rband_image[0].data)\n",
    "norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))\n",
    "ax.imshow(sdss_rband_image[0].data, cmap = 'gray_r', norm = norm, origin = 'lower')          \n",
    "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutout = Cutout2D(sdss_rband_image[0].data, pos, size, wcs=WCS(sdss_rband_image[0].header))\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1, projection=cutout.wcs)\n",
    "vmin,vmax = interval.get_limits(sdss_rband_image[0].data)\n",
    "norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))\n",
    "ax.imshow(cutout.data, cmap = 'gray_r', norm = norm, origin = 'lower')          \n",
    "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 15. Try looping over all positions and plotting multiwavelength cutouts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warning: this cell takes a long time to run! We limit it to the first three galaxies only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick the first 3 galaxies.\n",
    "galaxy_subset = galaxies[0:3]\n",
    "\n",
    "# For each galaxy,\n",
    "for galaxy in galaxy_subset:\n",
    "\n",
    "    # Establish the position.\n",
    "    ra = galaxy['RA']\n",
    "    dec = galaxy['DEC']\n",
    "    pos = SkyCoord(ra, dec, unit = 'deg')  \n",
    "    \n",
    "    # Set up the plot for this position.\n",
    "    fig = plt.figure(figsize=(20,6))\n",
    "    plt.suptitle('POSITION = ' + str(ra) + ', ' + str(dec), fontsize=16)\n",
    "\n",
    "    # GALEX\n",
    "    \n",
    "    # Find the GALEX images that overlap the position.\n",
    "    galex_image_table = galex_image_service.search(pos=pos, size=0.25)\n",
    "    \n",
    "    # Find the GALEX All-Sky Image Survey (AIS) Near-UV FITS coadd. \n",
    "    for i in range(len(galex_image_table)):\n",
    "        if (('image/fits' in galex_image_table[i].format) and\n",
    "            (galex_image_table['enrValue'][0]==2.35e-07)):\n",
    "            break\n",
    "    galex_image_record = galex_image_table[i]\n",
    " \n",
    "    # Create a cutout.\n",
    "    file_name = download_file(galex_image_record.getdataurl(), cache=True)  \n",
    "    gimage=fits.open(file_name)\n",
    "    galex_cutout = Cutout2D(gimage[0].data, pos, size, wcs=WCS(gimage[0].header))\n",
    "    \n",
    "    # Plot the cutout in the first position of a 1x3 (rowsxcols) grid.\n",
    "    ax = fig.add_subplot(1, 3, 1, projection=galex_cutout.wcs)\n",
    "    ax.set_title(galex_image_record.title)\n",
    "    ax.imshow(galex_cutout.data, cmap='gray_r', origin='lower', vmin = 0.0, vmax = 0.01)\n",
    "    ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n",
    "    \n",
    "    # SDSS\n",
    "    \n",
    "    # Find the SDSS images that overlap the position.\n",
    "    sdss_image_table = sdss_image_service.search(pos=pos, size=0, format='')\n",
    "    \n",
    "    # Find the first SDSS r-band image.\n",
    "    for sdss_rband_record in sdss_image_table:\n",
    "        if 'Sloan Digital Sky Survey - Filter r' in sdss_rband_record.title:\n",
    "            break\n",
    "    file_name = download_file(sdss_rband_record.getdataurl(), cache=True)  \n",
    "    sdss_rband_image=fits.open(file_name)\n",
    "    \n",
    "    # Create a cutout.\n",
    "    sdss_cutout = Cutout2D(sdss_rband_image[0].data, pos, size,\n",
    "                           wcs=WCS(sdss_rband_image[0].header))\n",
    "    \n",
    "    # Plot the cutout in the second position of a 1x3 grid.\n",
    "    vmin,vmax = interval.get_limits(sdss_cutout.data)\n",
    "    norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))\n",
    "    ax = fig.add_subplot(1, 3, 2, projection=sdss_cutout.wcs)\n",
    "    ax.imshow(sdss_cutout.data, cmap = 'gray_r', norm = norm, origin = 'lower')          \n",
    "    ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n",
    "    ax.set_title(sdss_rband_record.title)\n",
    "    \n",
    "    # AllWISE\n",
    "    \n",
    "    # Find the AllWISE images that overlap the position.\n",
    "    allwise_image_table = allwise_image_service.search(pos=pos, size=0)\n",
    "    \n",
    "    # Find the AllWISE W1 channel images.\n",
    "    for allwise_image_record in allwise_image_table:\n",
    "        if 'W1' in allwise_image_record.bandpass_id:\n",
    "            break\n",
    "    file_name = download_file(allwise_image_record.getdataurl(), cache=True)  \n",
    "    allwise_w1_image=fits.open(file_name)\n",
    "\n",
    "    \n",
    "    # Create a cutout.\n",
    "    allwise_cutout = Cutout2D(allwise_w1_image[0].data, pos, (size, size),\n",
    "                              wcs=WCS(allwise_w1_image[0].header))\n",
    "\n",
    "    # Plot the cutout in the third position of a 1x3 grid.\n",
    "    ax = fig.add_subplot(1, 3, 3, projection=allwise_cutout.wcs)\n",
    "    ax.imshow(allwise_cutout.data, cmap='gray_r', origin='lower', vmax = 10)\n",
    "    ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n",
    "    ax.set_title(allwise_image_record.title)\n"
   ]
  }
 ],
 "metadata": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}