HiPS_and_MOC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced usage of HiPS and MOCs\n",
    "## Exploring large catalogs within non-trivial spatial coverage, defined by brightness cuts and/or the availability of additional data sets\n",
    "Originally by \n",
    "\n",
    "**Caroline Bot and Thomas Boch, CDS**\n",
    "\n",
    "as published on http://www.euro-vo.org/?q=science/scientific-tutorials\n",
    "converted to jupyter notebook by \n",
    "\n",
    "**Katharina Lutz, Thomas Boch, Matthieu Baumann, CDS**\n",
    "\n",
    "This tutorial was first presented at \"Detecting the unexpected, Discovery in the Era of Astronomically Big Data\" Space Telsecope Science Instite, February 27 - March 2, 2017. It was updated to Aladin v10 in June 2017 and to Gaia DR2 in June 2018. \n",
    "\n",
    "This is a hands-on tutorial demonstrating an advanced useage of Hierarchical Progressive Suerveys (HiPS) and Multi-Order Coverage (MOC) maps in Aladin. Using this tutorial you will learn to handle a problem like: \n",
    ">\"I have a set of images, I would like to select regions in my observations that are above a given threshold in another survey (e.g. at low extinction), retrieve objects from very large catalogues (e.g. Gaia + WISE) in these non-trivial shapes and not-necessarily-connected regions and combine them to visualise some quantities (e.g. colour-colour diagram). \""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Before we get started\n",
    "The original tutorial was using Aladin (https://aladin.u-strasbg.fr/) and TOPCAT (http://www.star.bris.ac.uk/~mbt/topcat/), in this version of the tutorial we will use python packages and Aladin to accomplish the same tasks. The python packages include MOCpy (https://github.com/cds-astro/mocpy), astroquery (https://astroquery.readthedocs.io/en/latest/), pyVO (https://pyvo.readthedocs.io/en/latest/), ipyaladin (https://github.com/cds-astro/ipyaladin) as well as astropy and matplotlib. \n",
    "\n",
    "Now let's go!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings \n",
    "import glob\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "from astropy.coordinates import SkyCoord, Angle, match_coordinates_sky\n",
    "import astropy.wcs as astropy_wcs\n",
    "import astropy.units as u\n",
    "import mocpy\n",
    "import pyvo\n",
    "import healpy as hp\n",
    "import cdshealpix\n",
    "from astroquery.vizier import Vizier\n",
    "\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find all Short-Red images associated to the MASH catalog using the VizieR associated data service\n",
    "The VizieR service at CDS allows you to look for astronomical catalogues that have been published in the literature. Among these catalogues, some contain data associated to the publications and the tables therein. These data can be browsed and explored through the VizieR associated data service, which is linked to the traditional VizieR table service. In the current example we are looking for images associated with the MASH catalogue of planetary nebulae (Parker et al 2006-2008). The MASH fits files are cut-outs extracted from a larger H$\\alpha$ and Short Red survey and can be best described as a set of regions of interstes around planetary nebulae. \n",
    "\n",
    "To find VizieR associated data, we use the Table Access Protocol (TAP) with the VizieR endpoint. Through the VizieR TAP endpoint we can search for tables, content of tables and for information on associated data. \n",
    "\n",
    "First we search for the MASH catalogue:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tap_vizier = pyvo.dal.TAPService('http://tapvizier.u-strasbg.fr/TAPVizieR/tap')\n",
    "mash_catalogues = tap_vizier.search(\"SELECT  *  FROM tables \" + \n",
    "                                    \"WHERE description LIKE '%MASH%Parker%'\").to_table()\n",
    "mash_catalogues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this tutorial, we are interested in the tables belonging to the catalogues `V/127A`, this includes tables `V/127A/mash1` and `V/127A/mash2`. To have a look at the content of these tables we do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mash1_head = tap_vizier.search(\"SELECT TOP 5 * FROM \\\"V/127A/mash1\\\" \").to_table()\n",
    "mash1_head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the last column of this table is called `AssocData` and contains the entry `fits`. If you look at this table on the VizieR web interface, you can download the associated fits file. Within this notebook, we query the `obscor` database to get the URLs to the fits file. Using the `astropy.io.fits` module, we can then open the fits files from their URLs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mash_fits = tap_vizier.search(\"SELECT TOP 5  *  FROM obscore \" + \n",
    "                              \"WHERE obs_collection='V/127A'\").to_table()\n",
    "mash_fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the result from this query proivides us with information of the fits files, which are associated with the MASH catalogues. In particular the column `access_url` provides us with the location of the data. To get the first image we could do:\n",
    "\n",
    "`image = fits.open(mash_fits['access_url'][0])`\n",
    "\n",
    "and then work on the image, plot it or save it to our machine. However downloading all the data and working with all the data, takes quite some time. So, for the purpose of this tutorial, **we prepared a subsample of 335 of these Short Red images that will run in a timely manner** but we encourage you to repeat with the full Short Red sample later. **The subsample is available at \n",
    "http://astro.u-strasbg.fr/~bot/BochBot.tar.gz** and in the Data Folder of this repository. So if you run this tutorial on Binder, you do not need to download anything. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a MOC of the MASH images\n",
    "The multi-order coverage (MOC) map of a set of images is a useful represation of the sky coverage of these images. MOCs can represent arbitrary patches on the sky, which do not need to be connected. The union or intersection of two MOCs can be calculated quickly and with small computational effort. Catalogues can be filterd by MOCs. \n",
    "\n",
    "Here we can create the MOC of the MASH images with the `MOCPy` module. Please note that the majority of the following cell is concerned with removing keywords from the fits file headers, which would otherwise hamper the MOC creation (because the underlying `astropy.wcs.WCS` module is confused by these additional DEPRECATED header keywords):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a list of all fits files \n",
    "mash_file_list = glob.glob('Data/MASH_Sample/*.fits')\n",
    "# If this cell block was run before, there are too many files in the \n",
    "# folder, so before we do anything, we clean up.\n",
    "for file_name in mash_file_list:\n",
    "    if '_mod' in file_name:\n",
    "        os.remove(file_name)\n",
    "mash_file_list = glob.glob('Data/MASH_Sample/*.fits')\n",
    "# Modify the header so that the MOC creation works smoothly\n",
    "for file_name in mash_file_list:\n",
    "    print(file_name)\n",
    "    ima = fits.open(file_name)\n",
    "    removeable_keywords = ['CDELT1', 'CDELT2', 'PC001001', 'PC001002',\n",
    "                           'PC002001', 'PC002002', 'CROTA2', 'XPIXELSZ',\n",
    "                           'YPIXELSZ']\n",
    "    for keyword in removeable_keywords:\n",
    "        try:\n",
    "            del ima[0].header[keyword]\n",
    "        except KeyError:\n",
    "            continue\n",
    "    ima.writeto(file_name.split('.fits')[0] + '_mod.fits', overwrite=True)\n",
    "    ima.close()\n",
    "# Get a list of the updated files and create the MOC    \n",
    "mash_file_list = glob.glob('Data/MASH_Sample/*_mod.fits')\n",
    "moc_mash = mocpy.MOC.from_fits_images(mash_file_list, max_norder=15)\n",
    "moc_mash.write('Data/mash_moc.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('darkgrid')\n",
    "fig = plt.figure(figsize=(15, 15))\n",
    "with mocpy.WCS(fig, fov=2.5 * u.deg,\n",
    "               center=SkyCoord('21:49:00', '-02:30:00', \n",
    "                               unit=(u.hourangle, u.deg), frame='galactic'),\n",
    "               coordsys=\"galactic\", rotation=Angle(0, u.degree),\n",
    "               projection=\"AIT\") as wcs:\n",
    "    ax = fig.add_axes([0.17, 0.17, 0.77, 0.77], projection=wcs)\n",
    "    moc_mash.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color=mpl.cm.magma(0.4))\n",
    "    ax.set_xlabel('RA')\n",
    "    ax.set_ylabel('Dec')\n",
    "    lon, lat = ax.coords[0], ax.coords[1]\n",
    "    lon.set_major_formatter('hh:mm:ss')\n",
    "    lat.set_major_formatter('dd:mm')\n",
    "    lon.set_ticklabel(exclude_overlapping=True)\n",
    "    lat.set_ticklabel(exclude_overlapping=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This figure only shows a small region on the Sky, but you can see how the MOC has arbitrary shapes and not all regions are connected. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load an archival extinction map and create the MOC of the low extinction regions\n",
    "Different works (e.g. Schlegel et al. 1998, Schlafly  & Finkbeiner 2011, Green et al. 2015 or Dobashi et al. 2013) have created extinction maps of the sky and these maps are publicly available. Some of these maps are all-sky maps, other have higher resolutions, or come from different methods... Such maps are available in healpix format (among others) on the LAMBDA website or the CADE website. For the purpose of this tutorial, we will download the well known all-sky extinction map from Schlegel et al. from the LAMBDA website, and define as a MOC the low extinction area for which 0<E(B-V)<0.4. \n",
    "The map is available from here: https://lambda.gsfc.nasa.gov/data/foregrounds/SFD/lambda_sfd_ebv.fits and we save it to disc. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ext_map = fits.open('https://lambda.gsfc.nasa.gov/data/foregrounds/SFD/' + \n",
    "                    'lambda_sfd_ebv.fits')\n",
    "ext_map.writeto('Data/Schlegel_extinction_map.fits', overwrite=True)\n",
    "ext_map.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are only interested in regions with low extinction. So our aim is to get a MOC of all regions, where the extinction values from the Schlegel et al. map are in the range from 0 to 0.5mag. The extinction map we got from the NASA webpage is in the Healpix format, which is a very useful and efficient presentation of all-sky maps. The Healpix tesselation is also used by the MOCs. So to get the MOC from the extinction map, we do the following. First we change the projection of the map and set the order (i.e. resolution) of the map, so that the resulting MOC will have the right shape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the map\n",
    "ext_map = hp.read_map('Data/Schlegel_extinction_map.fits', nest=False)\n",
    "# define the coordinate transformation\n",
    "r = hp.Rotator(coord=['G','C'])\n",
    "# do the coordinate and representation transformation\n",
    "ext_map_equatorial_ring = r.rotate_map_pixel(ext_map)\n",
    "ext_map_equatorial = hp.reorder(ext_map_equatorial_ring, r2n=True)\n",
    "# set the resolution\n",
    "nside = 512\n",
    "norder = hp.nside2order(nside)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we declare which pixel we want to use, i.e. all pixels with low extinction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indexes_pixels_to_keep = np.where((ext_map_equatorial<0.5) & \n",
    "                                  (ext_map_equatorial>0.0))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we brought the extinction map in the right shape, the indexes of the pixels that we want to keep are the indexes of the Healpix cells we are interested in. Thus the MOC is created by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moc_low_ext = mocpy.MOC.from_healpix_cells(indexes_pixels_to_keep, \n",
    "                                           np.full((len(indexes_pixels_to_keep,)), \n",
    "                                                    norder))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moc_low_ext.write('Data/low_extinction_moc.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5:  Find out which regions are covered by the MASH short red images and in the low extinction regions defined above\n",
    "\n",
    "To find out which regions of the sky are covered by the MASH sample and which are at low extinction, we build the intersection of the two MOCs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moc_intersection = moc_low_ext.intersection(moc_mash)\n",
    "moc_intersection.sky_fraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can visualise the coverage of the two MOCs and their intersection. The grey area is where the extinction is low, the blue area is the MASH coverage and the red area is the MASH coverage in low extinction regions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(111, figsize=(10, 8))\n",
    "with mocpy.WCS(fig, \n",
    "               fov=200 * u.deg,\n",
    "               center=SkyCoord(200, -20, unit='deg', frame='icrs'),\n",
    "               coordsys=\"icrs\",\n",
    "               rotation=Angle(0, u.degree),\n",
    "               projection=\"AIT\") as wcs:\n",
    "    ax = fig.add_subplot(1, 1, 1, projection=wcs)\n",
    "    moc_low_ext.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color=\"grey\")\n",
    "    moc_mash.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color=\"dodgerblue\")\n",
    "    moc_intersection.fill(ax=ax, wcs=wcs, alpha=0.5, fill=True, color=\"crimson\")\n",
    "    ax.set_xlabel('RA')\n",
    "    ax.set_ylabel('Dec')\n",
    "    lon, lat = ax.coords[0], ax.coords[1]\n",
    "    lon.set_major_formatter('hh:mm:ss')\n",
    "    lat.set_major_formatter('dd:mm')\n",
    "    lon.set_ticklabel(exclude_overlapping=True)\n",
    "    lat.set_ticklabel(exclude_overlapping=True)\n",
    "    lon.set_ticks(spacing=2 * u.hourangle)\n",
    "plt.grid(color=\"black\", linestyle=\"dotted\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Query the WISE Catalogue by MOC\n",
    "Without the usage of MOC, querying for sources in the low extinction regions covered by the MASH subsample would be tedious or even impossible. Indeed, one would need to load the whole catalog and make selections which would not be possible given the size of some catalogues. Alternatively, one would need to query the catalogue field by field, which would take time and several queries. Instead, here, we will use the power of MOC files to query large catalogs directly in the covered regions only. We will use the coverage of the low extinction and MASH-covered region to query for sources from the Gaia and WISE surveys, in this highly non-continuous and non-trivial shape areas.\n",
    "\n",
    "First let's see which Gaia and WISE catalogues are available on VizieR. We could as above use the TAP endpoint of VizieR or we use as shown below the `Vizier` module in the `astroquery` package. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_list_wise = Vizier.find_catalogs('WISE')\n",
    "for k, v in catalog_list_wise.items():\n",
    "    print(k, ': ', v.description)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_list_gaia = Vizier.find_catalogs('Gaia DR2')\n",
    "for k, v in catalog_list_gaia.items():\n",
    "    print(k, ': ', v.description)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For WISE we will want to use `II/311 :  WISE All-Sky Data Release (Cutri+ 2012)` and for Gaia `I/345 :  Gaia DR2 (Gaia Collaboration, 2018)`. Before we actually query the two table we just get a few sources for each of the tables to understand, which columns are available ect. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_wise = Vizier.get_catalogs('II/311')\n",
    "print(test_wise)\n",
    "test_wise[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_gaia = Vizier.get_catalogs('I/345')\n",
    "print(test_gaia)\n",
    "test_gaia[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you will see below, we only need coordinates, WISE photometry in the W1 and W2 band, and Gaia photometry in the Gaia G band. So we'll query the tables `II/311/wise` for WISE and `I/345/gaia2` for Gaia DR2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wise = moc_intersection.query_vizier_table('II/311/wise', max_rows=100000)\n",
    "wise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia = moc_intersection.query_vizier_table('I/345/gaia2', max_rows=100000)\n",
    "gaia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7: Cross-match Gaia and WISE sources in all fields\n",
    "We now want to find sources in the selected region (observed in the MASH regions of interests and at low extinction) that are common to the WISE and Gaia catalogues. To do so, we will perform a cross-match of the Gaia and the WISE catalogues. Alternatively, we could also use the CDS XMatch service via the corresponding `astroquery` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wise_coord = SkyCoord(ra=wise['ra'], dec=wise['dec'], unit=u.deg)\n",
    "gaia_coord = SkyCoord(ra=gaia['ra_epoch2000'], \n",
    "                      dec=gaia['dec_epoch2000'], unit=u.deg)\n",
    "idx, d2d, d3d = match_coordinates_sky(wise_coord, gaia_coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_sep = 1.0 * u.arcsec\n",
    "sep_constraint = d2d < max_sep\n",
    "wise_matches = wise[sep_constraint]\n",
    "gaia_matches = gaia[idx[sep_constraint]]\n",
    "match_cat = wise_matches['JNAME', 'ra', 'dec', 'W1mag', 'W2mag']\n",
    "match_cat['Gmag'] = gaia_matches['phot_g_mean_mag']\n",
    "match_cat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 8: Build a colour-colour diagram\n",
    "\n",
    "We now use the data we got from the cross-match to get a WISE/Gaia colour-colour diagram for all the soruces in the sky regions that is covered by the MASH survey and at low extinction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('darkgrid')\n",
    "fig = plt.figure(figsize=(10.0, 8.0))\n",
    "ax = fig.add_axes([0.17, 0.17, 0.77, 0.77])\n",
    "ax.plot(match_cat['W1mag'] - match_cat['W2mag'], \n",
    "        match_cat['Gmag'] - match_cat['W1mag'], linestyle='', marker='.')\n",
    "ax.set_xlabel('W1-W2 [mag]')\n",
    "ax.set_ylabel('G-W1 [mag]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}