CS_Image_Access.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Searching for and retrieving images\n",
    "\n",
    "In this notebook, we show how to search for and retrieve images from VO services using the Registry and the __[Simple Image Access](http://www.ivoa.net/documents/SIA/)__ (SIA) protocol.\n",
    "\n",
    "1. [Finding SIA resources from the Registry](#finding) \n",
    "2. [Using SIA to retrieve an image](#sia) \n",
    "3. [Viewing the resulting image](#3.-Viewing-the-resulting-image)\n",
    "4. [Example with multiple data services](#4.-Example-of-data-available-through-multiple-services)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**\\*Note:**  for all of these notebooks, the results depend on real-time queries.  Sometimes there are problems, either because a given service has changed, is undergoing maintenance, or the internet connectivity is having problems, etc.  Always retry a couple of times, come back later and try again, and only then send us the problem report to investigate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline  \n",
    "\n",
    "import pyvo as vo\n",
    "\n",
    "from astropy.io import fits\n",
    "import astropy.coordinates as coord\n",
    "# For downloading files\n",
    "from astropy.utils.data import download_file\n",
    "\n",
    "import aplpy\n",
    "from IPython.display import Image as ipImage, display\n",
    "\n",
    "# There are a number of relatively unimportant warnings that 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": [
    "<a id=\"finding\"></a>\n",
    "\n",
    "# 1.  Finding SIA resources from the Registry\n",
    "\n",
    "First, how do we find out what  services are available?  These are listed in a registry at STScI (__[see here](http://www.ivoa.net/documents/RegTAP/)__).  Our Registry function gives a simple interface for how to search for services.  \n",
    "\n",
    "Let's search for services providing images in the ultraviolet bands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uv_services = vo.regsearch(servicetype='image',waveband='uv')\n",
    "uv_services.to_table()['ivoid','short_name','res_title']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This returns an astropy table containing information about the services available.  We can then specify the service we want by using the corresponding row.  We'll repeat the search with additional qualifiers to isolate the row we want (note that in the keyword search the \"%\" character can be used as a wild card):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uvot_services = vo.regsearch(servicetype='image',waveband='uv',keywords=['swift'])\n",
    "uvot_services.to_table()['ivoid','short_name','res_title']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows us that the data we are interested in comes from the HEASARC's SkyView service, but the point of these VO tools is that you don't need to know that ahead of time or indeed to care where it comes from."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"sia\"></a>\n",
    "\n",
    "# 2. Using SIA to retrieve an image:\n",
    "\n",
    "Now we look for images of our favorite source.  See __[the SIA definition](http://www.ivoa.net/documents/WD/SIA/sia-20040524.html)__ for usage.  In short, you can specify the central position and the size (degrees as one or two floats for the RA, DEC directions).  It is up to the service to determine how to provide this. Optionally, you can limit it to the format you want, e.g., \"image/fits\" or \"image/png\" etc.  \n",
    "\n",
    "What is returned to you is not the image itself but a list of images available and how to access them.  This is easiest shown by example:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = coord.SkyCoord.from_name(\"m51\")\n",
    "\n",
    "im_table = uvot_services[0].search(pos=coords,size=0.2,format='image/jpeg')\n",
    "im_table.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the fields you're interested in, e.g., the URLs of the images made by skyview.  Note that specifying as we did SwiftUVOT, we get a number of different images, e.g., UVOT U, V, B, W1, W2, etc.  For each survey, there are two URLs, first the FITS IMAGE and second the JPEG.  \n",
    "\n",
    "Note that different services will return different column names, but all will have a column giving the URL to access the image.  Though it has different column names in different services, it can always be accessed through the `getdataurl` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = im_table[0].getdataurl()\n",
    "print(url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Viewing the resulting image\n",
    "\n",
    "\n",
    "### JPG images\n",
    "\n",
    "Since we have asked for JPEG images, we can display an image in python easily by using its URL. Each row of the result has a getdataurl() method, and you can then hand the URL to an image displayer such as IPython.display:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = ipImage(url=im_table[0].getdataurl())\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fits files\n",
    "\n",
    "Or download the FITS image and display it with imshow, or aplpy.\n",
    "\n",
    "(This often errors off with a time out message.  Just try it again, possibly a couple of times.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Do the search again asking for FITS\n",
    "im_table = uvot_services[0].search(pos=coords,size=0.2,format='image/fits')\n",
    "\n",
    "#  Hand the url of the first result to fits.open()\n",
    "hdu_list = fits.open(im_table[0].getdataurl())\n",
    "hdu_list.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using imshow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(hdu_list[0].data, cmap='gray', origin='lower',vmax=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using aplpy\n",
    "\n",
    "This package knows how to read the header keywords in WCS to display the correct coordinate axes.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gc = aplpy.FITSFigure(hdu_list,figsize=(5, 5))\n",
    "gc.show_grayscale(stretch='log', vmax=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Example of data available through multiple services\n",
    "\n",
    "Suppose we want Sloan DSS data.  A generic query finds us a number of possibilities (note that this doesn't work for keywords=['sdss'];  be flexible and try several search terms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "services = vo.regsearch(servicetype='image', keywords=['sloan'], waveband='optical')\n",
    "services.to_table()[np.where(np.isin(services.to_table()['short_name'], b'SDSSDR7'))]['ivoid', 'short_name']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So one of these is served by SDSS's SkyServer and the other by HEASARC's SkyView.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using HEASARC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heasarc_dr7_service = [s for s in services if 'SDSSDR7' in s.short_name and 'heasarc' in s.ivoid][0]\n",
    "\n",
    "sdss_table_heasarc = heasarc_dr7_service.search(pos=coords,size=0.2,format='image/fits')\n",
    "sdss_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",
    "# hdu_list = sdss_table_heasarc[0].getdataobj()\n",
    "## But if you might run this notebook repeatedly with limited bandwidth, \n",
    "##  download it once and cache it.  \n",
    "\n",
    "#  Get the filter g version\n",
    "file_name=download_file(sdss_table_heasarc[0].getdataurl(),cache=True)\n",
    "hdu_list = fits.open(file_name)\n",
    "\n",
    "plt.imshow(hdu_list[0].data, cmap='gray', origin='lower', vmax=1200,vmin=1010)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using SDSS SkyServer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jhu_dr7_service = [s for s in services if 'SDSSDR7' in s.short_name and 'jhu' in s.ivoid][0]\n",
    "\n",
    "# Note: jhu_dr7_service access url has hard-wired \"format=image/fits\". \n",
    "# If you specify anythign else, it errors. If you specify nothing, \n",
    "# then the search() method puts \"format=all\", which errors. So specify empty string for now.\n",
    "sdss_table_jhu=jhu_dr7_service.search(pos=coords,size=0.2, format='')\n",
    "sdss_table_jhu.to_table().show_in_notebook(display_length = 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Get the filter g version\n",
    "file_name=download_file(sdss_table_jhu[1].getdataurl(),cache=True)\n",
    "hdu_list = fits.open(file_name)\n",
    "\n",
    "plt.imshow(hdu_list[0].data, cmap='gray', origin='lower',vmax=1200,vmin=1010)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that SkyView is just getting images by using the SIAP internally to get the data from the SDSS service.  The point of the VO protocols is that you don't need to know where the data are coming from.  But they can be processed differently.  "
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}