Handling gravitational-wave sky maps with Multi-Order Coverage (v2)

This document explains how gravitational wave sky maps can be easily and efficiently visualized and processed using Multi-Order Coverage (MOC) maps based on HEALPix sky tessellation. We compute the MOC region at a given probability level and subsequently, we query databases for retrieving objects whose position falls within this MOC map, we show operations between MOC maps and we check if a source is localized in a fixed level of probability.

For this tutorial we use the simulated sky maps from The First Two Years of Electromagnetic Follow-Up with Advanced LIGO and Virgo for compact binary Coalescence (CBC) sources. The sky maps for burst events are simulated from Localization of short duration gravitational-wave transients with the early advanced LIGO and Virgo detectors. A gallery of these sky maps is shown here.

We provide sample code in Python; you can download this document and run the code samples in IPython Notebook (now known as the Jupyter Notebook). The results are displayed in real time in Aladin Sky Atlas which is embedded in the document.

Table of Contents

  1. Installation and initialization
    1. Sky map visualization with Aladin
    2. Python packages
  2. Working with a gravitational-wave sky maps using MOC
    1. CBC sky maps
    2. Event id 18951
  3. MOC script commands
  4. Determining if an object falls into a specific level of probability
    1. Determining in which level of probability the source falls
  5. Query Catalogs from MOCs
    1. Query a single catalog
    2. Ranked list of galaxies in 3D sky map
    3. Queries running simultaneously
  6. Operation between MOC maps
    1. Intersection between MOC skymaps and VizieR tables footprints
  7. Utility
    1. Interactive MOCs

1. Installation and initialization

1.A Sky map visualization with Aladin

The interoperability between the Aladin Sky Atlas and the Python function outputs is obtained using SAMP (Simple Application Messaging Protocol). The best performance is achieved by installing Aladin Desktop and the Aladin java applet is embedded in the notebook. Anyhow, a detach button is supplied to disjoin the Aladin window.

As the code runs, the results are displayed in real time in the Aladin planes. Before running the notebook check the SAMP connection: Aladin bar --> Interop --> Connect with SAMP.

In [ ]:
from IPython.core.display import HTML
HTML('<iframe src=http://aladin.u-strasbg.fr/java/nph-aladin.pl width=770 height=930></iframe>')

1.B Python packages

We will need healpy for reading the probability sky map files and mocpy for parsing and manipulating MOCs.

For this tutorial, we will also use the astropy.vo astropy subpackage for accessing to the Virtual Observatory (VO) services. The astropy.utils subpackage contains general-purpose utility functions, astropy.table provides functionality for storing and manipulating heterogeneous tables, numpy for scientific computing with Python, urlparse is a standard interface to break Uniform Resource Locator (URL) strings up in components (the urlparse module is renamed to urllib.parse in Python 3) and future for running under Python 2 and Python 3. The json files are manipulated using json package. The last section makes use of ipywidgets; they are interactive HTML widgets for Jupyter notebooks and the IPython kernel.

SAMP is a protocol that is used by a number of other tools such as TOPCAT, SAO Ds9, and Aladin; it is possible to send and receive data to and from these tools. Here the Aladin script commands are converted in Python strings and sent via SAMPIntegratedClient class.

In [ ]:
from __future__ import print_function

from astropy.vo.samp import SAMPIntegratedClient
import os.path
import sys

class AladinViaSAMP(object):

    def __init__(self):
        self._client = SAMPIntegratedClient()
        
    def send_file(self, infile=str()):
        """Sending a file (image or table) to Aladin Sky Atlas using the SAMPIntegratedClient class.
             http://docs.astropy.org/en/stable/vo/samp/example_table_image.html
        """   
     
        self._client.connect()
        params = {}  
        
        if sys.version > '3':
            import urllib.parse 
            params[ "url" ] = urllib.parse.urljoin( 'file:', os.path.abspath( infile ) )
        else:
            import urlparse
            params[ "url" ] = urlparse.urljoin( 'file:', os.path.abspath( infile ) )
               
        message = {}
        message[ "samp.mtype" ] = "image.load.fits"
        message[ "samp.params" ] = params
     
        self._client.notify_all(message)
        self._client.disconnect()
        
    def send_script_command(self, script=str()):
        """Sending a script to Aladin Sky Atlas using the SAMPIntegratedClient class.
           http://docs.astropy.org/en/stable/vo/samp/example_table_image.html
         """

        self._client.connect()

        params = {}
        message = {} 
        message[ "samp.mtype" ] = "script.aladin.send"
        message[ "samp.params" ] = { "script" : script }  

        self._client.notify_all(message)
        self._client.disconnect()


class AladinScriptCommands(AladinViaSAMP):
    """A set of the main script commands for Aladin console.
        http://aladin.u-strasbg.fr/java/AladinScriptManual.gml"""

    def cview(self, url):
        """Creation of view: url."""
    
        cview_str = 'cview' + ' ' + url     
        return self.send_script_command(cview_str)

    def draw_circle(self, x, y, size = '10arcmin'):
        """Draw circle (x, y, size)."""
     
        position = [ x, y ]
        position = ' , '.join(map(str, position))  
        draw_circle_str = 'draw red circle' + '( ' + position + ', ' + size + ')'

        return self.send_script_command(draw_circle_str)

    def draw_line(self, line_values):
        """Draw line: .line(x1,y1,x2,y2,...[,text])."""
    
        draw_line_str = 'draw' + ' ' + 'line' + ' ' + line_values   
        return self.send_script_command(draw_line_str)

    def draw_newtool(self, name):
        """Create manually a new plane: draw newtool(name)."""
    
        draw_newtool_str = 'draw' + ' ' + 'newtool' + ' ' + name    
        return self.send_script_command(draw_newtool_str)

    def draw_string(self, x, y, text):
        """Draw string (x, y, text)."""

        position = [x, y]
        position = ' , '.join(map(str, position))
        draw_string_str = 'draw string' + '( ' + position + ', ' + text + ')'

        return self.send_script_command(draw_string_str)

    def draw_string_float(self, x, y, number):
        """Draw string (x, y, number)."""
     
        position = [ x, y ]
        position = ' , '.join(map(str, position))  
        draw_string_number = 'draw string' + '( ' + position +','+str(('% .1e' % number))+'%)'
     
        return send_script_command(self, draw_string_number)

    def get_FoV(self, x, y):
        """P_ra_dec = get FoV(pointing)."""
     
        position = [x, y] 
        position = '  '.join(map(str, position))

        plane_name = 'P:'+ str(x) + ',' + str(y) 
        get_fov_str =  plane_name + '= get FoV(pointing)' + ' ' + position

        return self.send_script_command (get_fov_str)

    def get_VizieR(self, catalog):
        """get VizieR(catalog,allsky)."""  
     
        get_vizier_str = 'get VizieR(' + catalog + ',' + 'allsky' + ')' 
        return self.send_script_command(get_vizier_str)

    def rename(self, plane):
        """Rename plane."""
        
        rename_plane_str = 'rename' + ' ' + plane      
        return self.send_script_command(rename_plane_str)

    def remove_FoV(self, x, y):
        """Remove Field of View"""

        position = [x, y]
        position = '  '.join(map(str, position))
        
        plane_name = 'P:'+ str(x) + ',' + str(y)
        remove_fov_str= 'rm' + ' ' + plane_name

        return self.send_script_command(remove_fov_str)

    def set_planeID(self, name):
        """Set planeID = todo"""

        set_plane_str = 'set' + ' ' + 'planeID=' + name
        return self.send_script_command(set_plane_str)

    def remove(self, plane):
        """Remove plane"""

        rm_str = 'rm' + ' ' + plane
        return self.send_script_command(rm_str)

    def cmoc(self, threshold, skymap, plane_name):
        """MOC extraction"""

        cmoc_str = plane_name+'='+'cmoc' + ' ' + '-threshold=' + str(threshold) + ' ' + skymap
          
        return self.send_script_command(cmoc_str)

    def set_moc(self, plane):
        """Set the properties of the MOC plane"""
        
        set_srt = 'set'+' '+plane +' '+'drawing=+perimeter,-border,-fill'

        return self.send_script_command(set_srt)

    def md(self, name):
        """Create a folder on the top of the Aladin stack"""

        md_srt = 'md'+ ' '+name

        return self.send_script_command(md_srt)

    def mv(self, plane, folder):
        """Move planes or views"""

        mv_srt = 'mv'+ ' ' + plane + ' ' + folder

        return self.send_script_command(mv_srt)

    def moc_inter(self, moc_1, moc_2, name):
        """MOC interpolation"""

        cmoc_inter_str = name + '=' + '' +'cmoc' +' ' + '-inter' + ' ' + moc_1 + ' ' + moc_2

        return self.send_script_command(cmoc_inter_str)

    def zoom(self, factor):
        """Change the zoom factor"""

        zoom_str = 'zoom' + ' ' + factor

        return self.send_script_command(zoom_str)

    def location(self, ra, dec):
        """Get sky position"""

        location_str = ra + ' ' + dec

        return self.send_script_command(location_str)  
    
    def get_hips (self, catalog):    
        """Call a remote image or tabular data server. """

        hips = 'get' + ' ' + 'hips(' + catalog + ')'
        
        return self.send_script_command(hips)

    def rm_all(self):
        """Remove all planes."""
    
        rm_all = 'rm -all'
    
        return self.send_script_command(rm_all)

The functions plot_interpolated_contours_from and _lvc_json are dedicated to manipulate .json files containing the contour plots of a probability sky map as shown in the Skymap Viewer. The interpolant is given in the healpy library using the four nearest points healpy.get_interp_val. The MOC represents an alternative approach to this. The two approaches are compared with one another in the Aladin folders. The first ones are referred as interpolated contour plots and the second ones as MOC contour plots or MOC confidence regions.

In [ ]:
import json
import sys
import os.path

aladin = AladinScriptCommands() # init.


def plot_interpolated_contours_from(json_link):
    """Plotting interpolated contour lines from an url or a local file.
    
    Input parameters
    ----------------
    json_link : string
              file/link of a LVC contour plot in json format.
    
    Return
    ------
    the probability levels are displayed in the Aladin Planes."""
    
    aladin.md("contour_plot") # creating a stack folder
    
    # download the json file from "url" and save it locally under "contour.json"
    if sys.version < '3':
        import urllib
        jsonfile = urllib.URLopener()
        jsonfile.retrieve( json_link, "contour.json" )
    
    else:
        import urllib.request
        urllib.request.urlretrieve( json_link, "contour.json" )
     
    with open( 'contour.json' ) as data_file:
        data = json.load( data_file )

    contour_pieces = len( data[ 'contours' ] )

    percentile = ('10-percentile','20-percentile','30-percentile','40-percentile',
                  '50-percentile','60-percentile','70-percentile',
                  '80-percentile','90-percentile')

    for percentile_json in percentile:
        aladin.draw_newtool ( percentile_json )
        _lvc_json( data, contour_pieces, percentile_json )
        aladin.mv(percentile_json, "contour_plot") #move in folder        
        
def _lvc_json(data, contour_pieces, percentile_json):
    """Managing the contour lines in a LVC json file."""
    
    i = 0
    for i in range(0, contour_pieces):
        contour = data[ 'contours' ][i]
        percentile = contour[ 'name' ]

        if percentile == percentile_json:
            values = contour[ 'coords' ]

            # sending to Aladin plane
            line = ( str( values ).replace('[' , '' ).replace(']' , '') )
            aladin.draw_line ( line )
In [ ]:
from __future__ import print_function

import healpy as hp
import numpy as np
from mocpy import MOC

from math import log

class MOC_confidence_region(object):
    """
    Multi-Order coverage map (MOC) of sky areas enclosed within a contour plot at a given confidence level.
    """

    def read_skymap(self, infile):
        """Reading healpix skymap.
        
        Input parameters
        ----------------
        infile : string
              LVC probability sky localization in healpix format
              
        Return
        -------
        hpx : list
            1D array of values (probability stored in each pixel)
        nside : int
           skymap resolution
        """      
        
        self.hpx = hp.read_map(infile, verbose = False)
        npix = len(self.hpx)
        self.nside = hp.npix2nside(npix)
        
        return self.hpx, self.nside
 
    def ipixs_in_percentage(self, percentage):
        """Finding ipix indices confined in a given percentage.
        
        Input parameters
        ----------------
        percentage : float
                 fractional percentage from 0 to 1  
        
        Return
        ------- 
        ipixs : list
              indices of pixels
        """
        
        sort = sorted(self.hpx, reverse = True)
        cumsum = np.cumsum(sort)
        index, value = min(enumerate(cumsum), key = lambda x: abs( x[1] - percentage ))

        index_hpx = range(0, len( self.hpx ))
        hpx_index = np.c_[self.hpx, index_hpx]

        sort_2array = sorted(hpx_index, key = lambda x: x[0], reverse = True)
        value_contour = sort_2array[0:index]

        j = 1 
        table_ipix_contour = [ ]

        for i in range (0, len(value_contour)):
            ipix_contour = int(value_contour[i][j])
            table_ipix_contour.append(ipix_contour)
    
        self.ipixs = table_ipix_contour
          
        return self.ipixs
     
    def sky_coords(self):
        """Creating an astropy.table with RA[deg] and DEC[deg] ipix positions
        
        Return
        ------- 
        contour_ipix : 
                    sky coords in degrees
        """
       
        # from index to polar coordinates
        theta, phi = hp.pix2ang(self.nside, self.ipixs)

        # converting these to right ascension and declination in degrees
        ra = np.rad2deg(phi)
        dec = np.rad2deg(0.5 * np.pi - theta)
        
        # creating an astropy.table with RA[deg] and DEC[deg]
        from astropy.table import Table
        self.contour_ipix = Table([ra, dec], names = ('RA[deg]', 'DEC[deg]'), 
                             meta = {'ipix': 'ipix table'})
     
        return self.contour_ipix
    
    def moc_order(self):
        """Setting MOC order.
        
        Return
        ------- 
        moc_order : int
              
        """       
        
        self.moc_order = int(log( self.nside, 2))
     
        return self.moc_order

    def create_moc(self):
        """Creating a MOC map from the contour_ipix table."""
        
        self.moc = MOC.from_table(self.contour_ipix, 'RA[deg]', 'DEC[deg]',
                                  self.moc_order)

        return self.moc

    def write_moc(self, percentage, short_name):
        """Writing MOC file in fits format.
        
        Input parameters
        ----------------
        percentage : float
                 fractional percentage from 0 to 1 converted into a string
        short_name : str
                 file output
        """
        
        return self.moc.write(short_name + '_MOC_' + str(percentage), format = 'fits')
    
    def contour_plot(self, infile, percentage, short_name=''):
        """Creating/Writing a MOC contour region at a fixed level of probability.
        
        Input parameters
        ---------------
        infile : string
              LVC probability sky localization in healpix format
        percentage : float
                 fractional percentage from 0 to 1
        """
        
        self.read_skymap(infile)
        self.ipixs_in_percentage(percentage)
        self.sky_coords()
        self.moc_order()
        self.create_moc()

        return self.write_moc(percentage, short_name)
    
    
class LocalizeSource(MOC_confidence_region):
    """The class is designed to localize an astrophysical source inside a probability skymap."""
        
    def __init__(self):
        self.aladin = AladinScriptCommands()
        
    def in_skymap(self, infile, ra, dec, percentage, label=' ', show = True):
        """Checking if an object falls in a given probability level.
        
        Input parameters
        ---------------
        infile : string
                LVC probability sky localization in healpix format
        percentage : float
                fractional percentage from 0 to 1 
        ra, dec : float
                sky coordinates in degrees
        label : string
                name of the object (optional, by default = '')
        show = True
                show the MOC confidence region and the object in the Aladin planes;
                otherwise no (optional; by default = True )
        """
        
        self.read_skymap(infile)
        ipixs=self.ipixs_in_percentage(percentage)
        
        theta = 0.5 * np.pi - np.deg2rad(dec)
        phi = np.deg2rad(ra)
        ipix = hp.ang2pix(self.nside, theta, phi)
        
        is_there = ipix in ipixs

        if is_there == True:
            print ("The sky coord", "ra="+str(ra)+"°,","dec="+str(dec)+"°", "(label:" + label+")", \
                    "lies within the", str(percentage*100)+'%', "c.l.")
        else:
            print ("The sky coord", "ra="+str(ra)+"°,","dec="+str(dec)+"°", "(label:" + label+")", \
                    "is outside the", str(percentage*100)+'%', "c.l.")

        if show == True:
            self.aladin.draw_newtool("sources")
            self.aladin.draw_string(ra, dec, label)
            self.aladin.draw_circle(ra, dec, size = '10arcmin')
    
    def pinpoint(self, infile, ra, dec, from_percentage=0.1, to_percentage=1,
                 resolution_percentage=0.1, label=' ', show=True):
            
        """Find in which confidence level the source falls.
        
        Input parameters
        ---------------
        infile : string
            LVC probability sky localization in healpix format
        from_percentage : float
            fractional percentage from 0 to 1
        to_percentage : float
            fractional percentage from 0 to 1
        resolution_percentage : float
            fractional percentage from 0 to 1
        ra, dec : float
            sky coordinates in degrees
        label : string
            name of the object (optional, by default = '')
        show = True
            show the MOC confidence region and the object in the Aladin planes;
        otherwise no (optional; by default = True )
        """
            
        self.read_skymap(infile)

        # finding ipix
        theta = 0.5 * np.pi - np.deg2rad(dec)
        phi = np.deg2rad(ra)
        ipix = hp.ang2pix(self.nside, theta, phi)

        find="n"
        while from_percentage < to_percentage or find=="y":
            ipixs = self.ipixs_in_percentage(from_percentage)
            is_there = ipix in ipixs

            if is_there != True:
                print ("It is not localized within the", str(from_percentage*100)+'%', "c.l.")
                from_percentage = from_percentage + resolution_percentage
            else:
                find="y"      
                print ("The sky coord", "ra="+str(ra)+"°,","dec="+str(dec)+"°", "(label:" + label+")" \
                       "lies within the", str(from_percentage*100)+'%', "c.l.")
                    
                self.aladin.draw_newtool("sources")
                self.aladin.draw_string(ra, dec, label)
                self.aladin.draw_circle(ra, dec, size = '10arcmin')
                    
                return find

2. Working with a gravitational-wave sky maps using MOC

MOC is a multiscale mapping based on HEALPix sky tessellation. It is essentially a simple way to map irregular and complex sky regions into hierarchically grouped predefined cells. Each MOC cell is defined by two numbers: the hierarchy level (HEALPIX ORDER) and the pixel index (HEALPIX NPIX).The NUNIQ scheme defines an algorithm for packing an (ORDER, NPIX) pair into a single integer for compactness:

$$uniq = 4\times 4^{(order)} + npix$$

Caution. By reducing a map in only a single confidence region, the probability distribution within that region is irreversible lost; see also Essick et al. (2015).

2.A CBC sky maps

The sky maps for compact binary Coalescence (CBC) are simulated from two analysis pipelines: the rapid pipeline - Bayestar - and the computationally intensive pipelines Lalinference_ MCMC or Lalinference_Nest.

The rapid Bayesian position reconstruction code that will produce accurate sky maps less than a minute after any BNS merger detection. The LALINFERENCE_MCMC (van der Sluys et al. 2008b; Raymond et al. 2009), LALINFERENCE_NEST (Veitch & Vecchio 2010), and LALINFERENCE_BAMBI (Graff et al. 2012, 2013) stochastic samplers were also used to follow up a subset of detected GW events. Though these analyses are significantly more computationally costly than BAYESTAR, taking hours to days, they can provide improved sky location estimates when the GW signal is very weak in one detector, and also yield not just sky localization but the full multidimensional probability distribution describing the parameters of a circularized compact binary merger.

2.B Event id 18951

This is the event 18951 - LIGO Hanford and Livingston joint detection. It simulates the response from a binary neutron star event at 75 Mpc. Here, we map the sky region in which the 90% of probability is enclosed. The MOC region is shown in white shadow over a discrete set of contour plots; each line encloses a given percentage of probability level from 10% to 90% in step of 10% (more about the contour plots).

In [5]:
# selecting an event id (2015);
# http://www.ligo.org/scientists/first2years/
# or  https://losc.ligo.org/s/skymapViewer/F2Y.html
event_id = '18951' 

# bayestar sky map
skymap_pipeline = 'bayestar'

# setting enclosed probability percentage 
prob_percentage = 0.5

# loading the simulated CBC event id (2015)
from astropy.utils.data import download_file

url_id = ('http://www.ligo.org/scientists/first2years/2015/compare/'
          + event_id +'/'+skymap_pipeline +'.fits.gz')

pipeline_event = download_file(url_id, cache = True, timeout = 300)

# sending to Aladin plane
aladin.send_file(pipeline_event)
aladin.rename (skymap_pipeline + event_id)

# MOC extraction: 
#   area enclosed within a specific contour plot at a given confidence level
moc = MOC_confidence_region() # init.

moc.contour_plot(infile = pipeline_event, 
                 percentage = prob_percentage, 
                 short_name = skymap_pipeline + event_id)

# sending it to Aladin plane
aladin.send_file(skymap_pipeline + event_id + 
               '_MOC_' + str(prob_percentage))

aladin.rename (skymap_pipeline + event_id + 
               '_MOC_' + str(prob_percentage))

# loading the MOC confidence region at 90%
MOC_file = MOC.from_file( skymap_pipeline + event_id 
                         + '_MOC_' + str(prob_percentage) )

# square degrees in a whole sphere
from math import pi
square_degrees_sphere = (360.0**2)/pi

# printing area
area_sq2 = round( ( MOC_file.sky_fraction * square_degrees_sphere ), 1 )
print ( str( int( prob_percentage*100 ) )+'%' + ' area = ', area_sq2, 'sq. deg' )

# loading DSS colored - sky background
aladin.get_hips( "P/DSS2/color" )

# plotting the interpoleted contour plots (Skymap Viewer plot)
plot_interpolated_contours_from( 'https://losc.ligo.org/s/skymapViewer/json/skymaps/F2Y/'+event_id+'.json' )
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 1.3 and will be removed in a future version. Use argument "overwrite" instead. [astropy.utils.decorators]
50% area =  159.4 sq. deg

Blink image between the contour plot and the MOC extraction of 90% area. See the completed and updated results in the section Sky map visualization with Aladin.

In [ ]:
# removing all planes
aladin.rm_all()

3. MOC script commands

The same results can be obtained using the MOC script commands provided in the Aladin Beta Version > v9.504.

In [6]:
# using the selected skymap in 2.B 
aladin.send_file(pipeline_event)
aladin.rename (skymap_pipeline + event_id)

# creating a stack folder       
aladin.md('MOC_contour') 

# creating moc contour plots from 10% to 90% in step of 10% 
for i in range(10 , 100, 10):  
    aladin.cmoc((i/100.0), skymap_pipeline + event_id, 
                'moc'+str(i/100.0)+skymap_pipeline + event_id)

for i in range(10 , 100, 10): #move in folder
    aladin.mv('moc'+str(i/100.0) + skymap_pipeline + event_id,
              'MOC_contour') # moving in folder
    
for i in range(10 , 100, 10): # drawing method: perimeter
    aladin.set_moc('moc'+str(i/100.0)+skymap_pipeline + event_id) 
    
# loading DSS colored - sky background
aladin.get_hips( "P/DSS2/color" )

Generation of MOC contour plots from 0.1 to 0.9 probability threshold (enclosed probability) in step of 0.1 using Aladin Beta Version v9.504. See the partial results in the section Sky map visualization with Aladin.

4. Determining if an object falls into a specific level of probability

Is an object located within a given level of probability? To answer this question, we implemented the method ._in_skymap(infile, ra, dec, percentage, label=' ', show = True) in LocalizeSource class. The method takes as an input a probability skymap in healpix format (infile), a fraction of the enclosed probability (percentage), the sky coord of the target in degrees (ra, dec) with an appropriate tag (label, by default is an empty string). If Show = True, the target is displayed in the Aladin plane.

As an example, we consider a list of transients and define if they lie within the 90%. The sky positions are shown in the Aladin planes.

In [7]:
# setting input values
infile = pipeline_event # using the selected skymap in 2.B 
ra_transients = [346.83421, 336.46342, 328.25277] # deg unit
dec_transients = [+43.99271, +24.56405, +19.44063] # deg unit
label_transients = "Source_1", "Source_2", "Source_3"
percentage = 0.9

localize = LocalizeSource() # init.
for ra_transient, dec_transient, label_transient in zip(ra_transients, dec_transients, label_transients):
    localize.in_skymap(infile, ra_transient, dec_transient, percentage, label_transient, show = True)
The sky coord ra=346.83421°, dec=43.99271° (label:Source_1) is outside the 90.0% c.l.
The sky coord ra=336.46342°, dec=24.56405° (label:Source_2) is outside the 90.0% c.l.
The sky coord ra=328.25277°, dec=19.44063° (label:Source_3) lies within the 90.0% c.l.

The method ._in_skymap(infile, percentage, ra, dec, label=' ', show = True) determines if an object falls into a specific level of probability; by default the result is shown in the Aladin stack. See the completed and updated results in the section Sky map visualization with Aladin.

4.A Determining in which level of probability the source falls

The method pinpoint(self, infile, ra, dec, from_percentage=0.1, to_percentage=1, resolution_percentage=0.1, label=' ', show=True) in LocalizeSource class is designed to determine in which level of probability a source is localized. The method takes as an input a probability skymap in healpix format (infile), a range of the enclosed probability (from_percentage=0.1, to_percentage=1, resolution_percentage=0.1), the sky coord of the target in degrees (ra, dec) with an appropriate tag (label, by default is an empty string). If Show = True, the target is displayed in the Aladin plane. The localization accuracy is improved with increasing resolution_percentage.

We evaluate the position of Source_3 The sky coord ra=328.25277°, dec=19.44063° (label:Source_3) is inside the 90.0% MOC contour plot

In [8]:
localize = LocalizeSource() # init.

localize.pinpoint(infile=pipeline_event, ra=328.25277, dec=19.44063, from_percentage=0.2, to_percentage=1,
                  resolution_percentage=0.01, label='Source_3', show=True)
It is not localized within the 20.0% c.l.
It is not localized within the 21.0% c.l.
It is not localized within the 22.0% c.l.
It is not localized within the 23.0% c.l.
It is not localized within the 24.0% c.l.
The sky coord ra=328.25277°, dec=19.44063° (label:Source_3)lies within the 25.0% c.l.
Out[8]:
'y'
In [9]:
# removing all planes
aladin.rm_all()

5. Query Catalogs from MOCs

In this section, we show how MOC maps can be use to query catalog objects that falls into the sky map region. The MOCs of all VizieR tables footprints are available on line (about 16.000 tables) and can be queried simultaneously in few seconds.

5.A Query a single catalog

As an example, we query the Gravitational Wave Galaxy Catalog from the MOC map obtained in section 2.B Event id 18951. The source positions are displayed in red inside the MOC region drawn in white.

In [10]:
catalog = 'VII/267/gwgc' # selecting catalog
catalog_renamed = catalog.replace('/', '_')

# loading MOC confidence region obtained in 2.B 
from mocpy import MOC
moc = MOC.from_file( 'bayestar18951_MOC_0.9' )

# sending to Aladin plane
aladin.send_file( 'bayestar18951_MOC_0.9')
aladin.rename ( 'bayestar18951_MOC_0.9' )

# querying from MOC ignoring astropy.io.votable.exceptions
import warnings
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    table = moc.query_vizier_table( catalog, max_rows = 100000 ) 

# file output: votable format
table.write( catalog_renamed + 'MOC_query', format = 'votable', overwrite = True )

# sending to the Aladin plane 
aladin.send_file( catalog_renamed + 'MOC_query' )

# loading DSS colored for sky background
aladin.get_hips( "P/DSS2/color" )

Catalog query from a MOC probability map (90% conf. level). See the completed and updated results in the section Sky map visualization with Aladin

In [11]:
# removing all planes
aladin.rm_all()

5.B Ranked list of galaxies in 3D sky map

Singer et al. (2016) discuss a rapid algorithm for obtaining joint three-dimensional estimates of sky location and luminosity distance from observations of binary neutron star mergers with Advanced LIGO and Virgo.
They argued that combining the reconstructed volumes with positions and redshifts of possible host galaxies can provide a manageable list of targets to search for optical or infrared emission. The 2MASS Redshift Survey (2MRS) (Huchra et al. 2012) is downloaded for this purpose. In order to reduce the query time, the 2MRS is directly queried from the MOC region e.g in which the 90% of probability is enclosed. Here, the MOC map obtained in section 2.B Event id 18951 is used.

See here for more details on the sample python code provided in the online supplement to the Letter GOING THE DISTANCE: MAPPING HOST GALAXIES OF LIGO AND VIRGO SOURCES IN THREE DIMENSIONS USING LOCAL COSMOGRAPHY AND TARGETED FOLLOW-UP.

In [20]:
# downloading 3D HEALPix sky map
from astropy.utils.data import download_file
url = ('http://asd.gsfc.nasa.gov/Leo.Singer/'+
       'going-the-distance/2015/compare/18951/'+'bayestar.fits.gz')

filename = download_file(url, cache=True)

# reading HEALPix layers
import healpy as hp
prob, distmu, distsigma, distnorm = hp.read_map(filename, 
                                                field=[0, 1, 2, 3], verbose=False)

# HEALPix resolution 
npix = len(prob)
nside = hp.npix2nside(npix)

pixarea = hp.nside2pixarea(nside)

# Ranking list of galaxies from a MOC region
from mocpy import MOC
moc = MOC.from_file( 'bayestar18951_MOC_0.9' )

# sending to Aladin plane
aladin.send_file( 'bayestar18951_MOC_0.9')
aladin.rename ( 'bayestar18951_MOC_0.9' )

catalog = 'J/ApJS/199/26/table3'  # 2MASS Redshift Survey 

# querying from MOC ignoring astropy.io.votable.exceptions
import warnings
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    cat = moc.query_vizier_table( catalog, max_rows = 100000 )

import numpy as np
from scipy.special import gammaincinv

completeness = 0.5
alpha = -1.0
MK_star = -23.55

MK_max = MK_star + 2.5*np.log10(gammaincinv(alpha + 2, completeness))

# selecting only galaxies with positive redshifts and absolute
# magnitudes greater than M(max)
from astropy.cosmology import WMAP9 as cosmo

import astropy.units as u
import astropy.constants as c

z = (u.Quantity(cat['cz']) / c.c).to(u.dimensionless_unscaled)

MK = cat['Ktmag'] - cosmo.distmod(z)
keep = (z > 0) & (MK < MK_max)

cat = cat[keep]
z = z[keep]

# luminosity distance and HEALPix index of each galaxy
r = cosmo.luminosity_distance(z).to('Mpc').value

theta = 0.5*np.pi - cat['_DEJ2000'].to('rad').value
phi = cat['_RAJ2000'].to('rad').value
ipix = hp.ang2pix(nside, theta, phi)

# probability density per unit volume at the position of each galaxy
from scipy.stats import norm
dp_dV = prob[ipix]*distnorm[ipix]*norm(distmu[ipix], distsigma[ipix]).pdf(r) / pixarea

#sorting the galaxies by descending probability density
galaxies_in_moc = cat[np.flipud(np.argsort(dp_dV))][:]

# adding probability galaxy position to the catalog
from astropy.table import Column

dp_dV_sort = np.flipud(np.argsort(dp_dV))[:]
dp_dV_value = dp_dV[dp_dV_sort]

# rounding
dp_dV_value_round = []
dp_dV_value_round = ['{:.3e}'.format(i) for i in dp_dV_value]

probability_galaxy_position = Column(dp_dV_value_round, name = 'dp_dV')

galaxies_in_moc.add_column(probability_galaxy_position, index=0)
print (galaxies_in_moc['_RAJ2000', '_DEJ2000', 'Ktmag','dp_dV'])

# sending to Aladin plane the weighted catalog 
galaxies_in_moc.write( 'ranked_list_galaxies', format = 'votable', overwrite = True )
aladin.send_file( 'ranked_list_galaxies' )

# loading DSS colored - sky background
aladin.get_hips( "P/DSS2/color" )
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:54: RuntimeWarning: invalid value encountered in less
 _RAJ2000  _DEJ2000 Ktmag    dp_dV   
   deg       deg     mag             
--------- --------- ------ ----------
344.01190  36.36136  8.772  5.476e-05
343.81122  36.67177  9.958  5.167e-05
137.19089 -38.60788  9.566  4.673e-05
334.86545  29.39581  9.835  4.640e-05
359.81589  46.88923  9.307  4.150e-05
  0.00695  47.27456  9.499  4.060e-05
343.35303  37.17527 10.151  4.052e-05
137.01408 -37.56026 10.388  3.510e-05
353.71420  44.30398  9.733  3.277e-05
139.60875 -39.28682  9.968  3.236e-05
      ...       ...    ...        ...
191.14339 -69.19820 11.434  2.787e-58
328.78137  21.45712 11.280  6.323e-59
317.35141  -1.96921 11.473  9.233e-62
195.17738 -69.79794 11.357  1.760e-73
204.73236 -70.26720 10.858  5.169e-74
313.54910  -2.76770 11.686  3.482e-78
124.54511 -12.71006 11.409  2.059e-79
128.63509 -26.21134 11.189  5.328e-87
 10.31410  47.93468 10.919 2.518e-114
 38.51760  54.57246 11.472 6.155e-143
Length = 332 rows

2MASS Redshift Survey catalog; in the first column the probability density per unit volume at the position of each galaxy is reported. See the completed and updated results in the section Sky map visualization with Aladin

In [13]:
# removing all planes
aladin.rm_all()

5.C Queries running simultaneously

Here we simultaneously query 6 VizieR tables footprints. The following catalogs are selected: a) Gravitational Wave Galaxy Catalogue (White 2011); b) Compact Binary Coalescence Galaxy Catalog (Kopparapu et al., 2008); c) Catalogue of Rich Clusters of Galaxies (Abell et al., 1989); d) Northern Cluster Catalog (Gal et al., 2009); e) MCXC Meta-Catalogue X-ray galaxy Clusters (Piffaretti et al., 2011); f) 2MASS Redshift Survey (Huchra et al. 2012) Here the MOC map obtained in section 2.B Event id 18951 is used.

In [14]:
# selecting catalogs
catalogs = ['VII/267/gwgc','J/ApJ/675/1459/table1','VII/110A',
            'J/AJ/137/2981','J/A+A/534/A109','J/ApJS/199/26/table3'] 

# loading MOC confidence region obtained in 2.B
from mocpy import MOC
moc = MOC.from_file( 'bayestar18951_MOC_0.9') 

# sending to Aladin plane
aladin.send_file( 'bayestar18951_MOC_0.9')
aladin.rename ( 'bayestar18951_MOC_0.9' )

# querying from MOC ignoring astropy.io.votable.exceptions
import warnings
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    for catalog in catalogs:
        print ()
        table = moc.query_vizier_table( catalog, max_rows = 100000 ) 
        print (table)
        catalog_renamed = catalog.replace('/', '_')
        table.write( catalog_renamed + 'MOC_query', format = 'votable', overwrite = True )
        aladin.send_file( catalog_renamed + 'MOC_query' )
        
# loading DSS colored for sky background
aladin.get_hips( "P/DSS2/color" )
_RAJ2000 _DEJ2000   PGC        Name      RAJ2000  ... e_Dist Simbad NED LEDA
  deg      deg                             deg    ...  Mpc                  
-------- -------- ------- -------------- -------- ... ------ ------ --- ----
  0.0070  47.2745       2       UGC12889   0.0070 ...  10.87 Simbad NED LEDA
  0.0540  46.9651      18      PGC000018   0.0540 ...  11.60 Simbad NED LEDA
  1.1110  47.4903     304      PGC000304   1.1110 ...  11.39 Simbad NED LEDA
  1.3665  46.5436     380      PGC000380   1.3665 ...  10.79 Simbad NED LEDA
  1.6535  47.8789     499       UGC00048   1.6534 ...   9.42 Simbad NED LEDA
  1.7451  46.6379     540       UGC00055   1.7451 ...  11.12 Simbad NED LEDA
  1.8490  47.0407     574       UGC00061   1.8490 ...  11.40 Simbad NED LEDA
  1.8523  46.9888     575      PGC000575   1.8523 ...  11.03 Simbad NED LEDA
  2.0605  49.7685     623       UGC00070   2.0605 ...  16.54 Simbad NED LEDA
  2.3694  47.3558     676       UGC00085   2.3694 ...  11.15 Simbad NED LEDA
     ...      ...     ...            ...      ... ...    ...    ... ...  ...
203.0311 -69.4121 4614889 HIPASSJ1332-69 203.0311 ...   8.70 Simbad NED LEDA
136.6347 -35.4732 4614928 HIPASSJ0906-35 136.6347 ...  13.89 Simbad NED LEDA
132.2641 -33.0483 4614931 HIPASSJ0848-33 132.2641 ...   3.43 Simbad NED LEDA
128.5590 -18.7043 4614935 HIPASSJ0834-18 128.5590 ...  12.28 Simbad NED LEDA
136.4961 -36.0604 4663310 HIPASSJ0905-36 136.4961 ...   2.09 Simbad NED LEDA
136.7370 -34.3194 4668280 HIPASSJ0906-34 136.7370 ...   3.71 Simbad NED LEDA
135.2350 -37.1229 4668309 HIPASSJ0900-37 135.2350 ...   2.74 Simbad NED LEDA
136.9908 -37.2214      --          Pyxis 136.9908 ...     -- Simbad NED LEDA
194.8950 -70.8747      --        NGC4833 194.8950 ...     -- Simbad NED LEDA
322.4925  12.1669 2802701        NGC7078 322.4925 ...     -- Simbad NED LEDA
Length = 384 rows

_RAJ2000 _DEJ2000           Name          RAJ2000  ...    Lum     errB e_Dist
  deg      deg                            "h:m:s"  ... 1e+10 Sun  mag   kpc  
-------- -------- ----------------------- -------- ... --------- ----- ------
  9.7500  48.3333                 NGC0185 00 39 00 ...     0.006  0.42   0.10
  8.2750  48.5167                 NGC0147 00 33 06 ...     0.006  0.42   0.10
337.1075  30.2913                 NGC7292 22 28 26 ...     0.066  0.30   0.20
132.2750 -26.3215              ESO496-010 08 49 06 ...     0.022  0.42   0.25
351.1750  41.3333                UGC12588 23 24 42 ...     0.077  0.42   0.22
352.5000  41.0000                UGC12632 23 30 00 ...     0.147  0.42   0.22
146.3525 -48.1413              PGC2791980 09 45 25 ...     0.016  0.42   0.25
129.0625 -26.4093              ESO495-021 08 36 15 ...     0.229  0.42   0.25
350.5275  40.8455                 NGC7640 23 22 07 ...     1.259  0.30   0.20
156.1850 -54.7973              ESO168-002 10 24 44 ...     0.520  0.42   0.25
     ...      ...                     ...      ... ...       ...   ...    ...
340.3400  34.6208              PGC2052000 22 41 22 ...     0.056  0.42   0.25
339.4175  34.8462                UGC12121 22 37 40 ...     3.342  0.42   0.25
335.6500  28.4688               PGC068685 22 22 36 ...     2.780  0.42   0.25
315.2825   0.1953                 NGC7001 21 01 08 ...     8.872  0.42   0.25
209.1550 -70.9275              ESO066-003 13 56 37 ...     9.638  0.42   0.25
341.2675  33.9962                UGC12179 22 45 04 ...     5.152  0.42   0.25
315.2400   0.3173 SDSSJ210057.49-001902.7 21 00 58 ...     0.337  0.42   0.25
314.8700   0.3547 SDSSJ205928.87+002116.8 20 59 29 ...     0.089  0.42   0.25
314.4825   0.1198              PGC1157607 20 57 56 ...     0.912  0.42   0.25
 20.1300  50.1445               PGC004829 01 20 31 ...     5.200  0.42   0.25
Length = 274 rows

_RAJ2000 _DEJ2000 ACO  RAB1950 DEB1950 ... Rich Dclass m10   _RA.icrs   _DE.icrs
  deg      deg         "h:m:s" "d:m:s" ...             mag   "h:m:s"    "d:m:s" 
-------- -------- ---- ------- ------- ... ---- ------ ---- ---------- ---------
  8.9646  49.6919   63 00 33.1  +49 25 ...    0      4 16.1 00 35 51.5 +49 41 31
324.0237  14.4411 2359 21 33.7  +14 13 ...    0      6 17.6 21 36 05.7 +14 26 28
328.3942  17.6697 2390 21 51.2  +17 26 ...    1      6 17.6 21 53 34.6 +17 40 11
330.2377  20.9576 2409 21 58.6  +20 43 ...    2      5 16.8 22 00 57.0 +20 57 27
346.6847  40.6539 2535 23 04.4  +40 23 ...    1      5 16.9 23 06 44.3 +40 39 14

 _RAJ2000  _DEJ2000 n_NSC      NSC       ...   Lopt    e_Lopt  __off_   beta 
   deg       deg                         ... [solLum] [solLum]  Mpc          
--------- --------- ----- -------------- ... -------- -------- ------ -------
318.27631   5.93010     X J211306+055548 ...    1.068    0.313   0.03  -23.20
318.58645   6.73020     S J211420+064348 ...    0.231    0.240   0.16   11.70
318.68603   3.48452     X J211444+032904 ...    0.753    0.226   0.26    2.00
318.68719   2.56750     X J211444+023403 ...    0.844    0.249   0.09   -0.70
318.80031   3.10298       J211512+030610 ...    0.140    0.157   0.31      --
318.99674   7.76171     S J211559+074542 ...    0.303    0.146   0.04    5.10
318.99794   7.45468     S J211559+072716 ...    0.333    0.169   0.10  -24.20
319.08391   4.35612     S J211620+042122 ...    0.262    0.239   0.11   20.00
319.21518   7.73074     S J211651+074350 ...    0.262    0.216   0.08   15.30
319.38867   6.88704     S J211733+065313 ...    0.437    0.186   0.25  -16.00
      ...       ...   ...            ... ...      ...      ...    ...     ...
324.18440   9.23196       J213644+091355 ...    0.282    0.226   0.17      --
324.21832  12.79893       J213652+124756 ...    0.000    0.000   0.00      --
324.24634  10.03088       J213659+100151 ...    0.286    0.331   0.01      --
324.43304  10.40843     S J213743+102430 ...    0.270    0.178   0.10   -7.10
324.50364  10.13910     S J213800+100820 ...    0.306    0.184   0.03  -41.30
324.87240  10.48039       J213929+102849 ...    0.361    0.256   0.12      --
324.98904  12.77389     S J213957+124626 ...    0.253    0.157   0.18   36.50
328.24368  17.67592       J215258+174033 ...    0.761    0.342   0.13      --
328.44750  17.64328     X J215347+173835 ...    0.425    0.155   0.04  -16.30
328.48599  17.25697     S J215356+171525 ...    0.223    0.167   0.59    4.60
Length = 45 rows

_RAJ2000 _DEJ2000     MCXC          OName       ...   R500  Notes Simbad
  deg      deg                                  ...   Mpc               
-------- -------- ------------ ---------------- ... ------- ----- ------
 16.9492  54.1419 J0107.7+5408 RXC J0107.7+5408 ...  1.0415       Simbad
126.6925 -20.1328 J0826.7-2007 RXC J0826.7-2007 ...  0.8460       Simbad
318.4746   2.5556 J2113.8+0233 RXC J2113.8+0233 ...  0.8197       Simbad
328.3979  17.6867 J2153.5+1741 RXC J2153.5+1741 ...  1.3554       Simbad
330.2200  20.9708 J2200.8+2058 RXC J2200.8+2058 ...  1.1217       Simbad
349.6637  42.9658 J2318.6+4257 RXC J2318.6+4257 ...  0.5969       Simbad
350.0600  41.7786 J2320.2+4146 RXC J2320.2+4146 ...  0.9592       Simbad

 _RAJ2000  _DEJ2000        ID         A  ... _2MX        SimbadName       NED
   deg       deg                         ...                                 
--------- --------- ---------------- --- ... ---- ----------------------- ---
339.26709  34.41592 22370410+3424573     ...  2MX 2MASX J22370410+3424573 NED
128.34517 -22.97367 08332284-2258252     ...  2MX 2MASX J08332284-2258252 NED
  9.74154  48.33738 00385796+4820145     ...  2MX 2MASX J00385796+4820145 NED
  8.30050  48.50874 00331212+4830314     ...  2MX 2MASX J00331212+4830314 NED
139.65295 -38.01009 09183669-3800363     ...  2MX 2MASX J09183669-3800363 NED
133.63544 -32.93748 08543250-3256149   z ...  2MX 2MASX J08543250-3256149 NED
126.04198 -18.77556 08241008-1846318     ...  2MX 2MASX J08241008-1846318 NED
121.40851 -11.42704 08053803-1125372     ...  2MX 2MASX J08053803-1125372 NED
344.01190  36.36136 22560284+3621411     ...  2MX 2MASX J22560284+3621411 NED
129.14565 -20.47047 08363496-2028136     ...  2MX 2MASX J08363496-2028136 NED
      ...       ...              ... ... ...  ...                     ... ...
145.55591 -45.37045 09421341-4522138     ...  2MX 2MASX J09421341-4522138 NED
137.94478 -40.38942 09114676-4023220     ...  2MX 2MASX J09114676-4023220 NED
179.61818 -67.57728 11582838-6734384     ...  2MX 2MASX J11582838-6734384 NED
187.15013 -69.51701 12283604-6931007     ...  2MX 2MASX J12283604-6931007 NED
131.77176 -27.79185 08470522-2747306     ...  2MX 2MASX J08470522-2747306 NED
136.70439 -34.16213 09064903-3409436     ...  2MX 2MASX J09064903-3409436 NED
327.88501  21.19614 21513240+2111461     ...  2MX 2MASX J21513240+2111461 NED
138.09196 -38.25584 09122205-3815211     ...  2MX 2MASX J09122205-3815211 NED
195.79318 -70.63913 13031040-7038212     ...  2MX 2MASX J13031040-7038212 NED
338.51093  33.81857 22340262+3349070     ...  2MX 2MASX J22340262+3349070 NED
Length = 506 rows

Multi-catalog query from a MOC probability map (90% conf. level). See the completed and updated results in the section Sky map visualization with Aladin

In [15]:
# removing all planes
aladin.rm_all()

6. Operation between MOC maps

6.A Intersection between MOC skymaps and VizieR tables footprints

The operations between the MOC maps (union, intersection, subtraction, difference) are extremely simple and fast (generally a few milliseconds) even for very complex sky regions. Here the intersection between the DSS coverage and the MOC sky map in 2.B Event id 18951 is shown.

In [16]:
# loading MOC confidence region obtained in 2.B
from mocpy import MOC
moc_1 = MOC.from_file( 'bayestar18951_MOC_0.9' ) 

# sending to Aladin plane
aladin.send_file( 'bayestar18951_MOC_0.9')
aladin.rename ( 'bayestar18951_MOC_0.9' )

# loading the MOC coverage map of SDSS Photometric Catalog (9)
from astropy.utils.data import download_file
url_id = 'http://alasky.u-strasbg.fr/footprints/tables/vizier/V_139_sdss9/MOC'
sdss9 = download_file( url_id, cache = True, timeout = 300 )
aladin.send_file( sdss9 )
aladin.rename ( 'sdss9_MOC' )

#load sdss9 MOC coverage
moc_2 = MOC.from_file( sdss9 ) 

# Intersection operation and writing file
inter = moc_1.intersection( moc_2 )
inter.write( 'intersection', format = 'fits')

#sending to Aladin plane
aladin.send_file( 'intersection' )
aladin.rename ( 'intersection' )

# loading DSS colored - sky background
aladin.get_hips( "P/DSS2/color" )

Intersection (in red) between the MOC map of SDSS Photometric Catalog (in blue) and MOC sky map at the 90% probability of event id 18951 (in white). See the completed and updated results in the section Sky map visualization with Aladin

In [17]:
# removing all planes
aladin.rm_all()

7. Utility

7.A Interactive MOCs

The probability MOC regions are automatically generated by adjusting the probability threshold slider. The 2.B Event id 18951 sky map is choosen in the example.

In [19]:
import ipywidgets as widgets
from ipywidgets import interact, fixed

# selecting an event id (2015);
# http://www.ligo.org/scientists/first2years/
event_id = '11762'

# bayestar sky map
skymap_pipeline = 'bayestar'

# loading the simulated CBC event id (2015)
from astropy.utils.data import download_file

url_id = 'http://www.ligo.org/scientists/first2years/2015/compare/'+event_id+'/'+skymap_pipeline+'.fits.gz'
pipeline_event = download_file( url_id, cache = True, timeout = 300 )

# sending to Aladin plane
aladin.send_file ( pipeline_event )
aladin.rename ( skymap_pipeline + event_id )

#slider MOC production
interact( MOC_confidence_region, infile = pipeline_event, 
         percentage = (0.1, 0.9, 0.1), short_name = fixed( skymap_pipeline + event_id ) )

# loading DSS colored for sky background
aladin.get_hips( "P/DSS2/color" )

# plotting the interpoleted contour plots (Skymap Viewer plot)
plot_interpolated_contours_from( 'https://losc.ligo.org/s/skymapViewer/json/skymaps/F2Y/'+event_id+'.json' )



# basic MOC function for slider
def MOC_confidence_region(infile, percentage, short_name = ' '):   
    """  Multi-Order coverage map (MOC) of sky area enclosed within a contour plot at a given confidence level.
    
    Input:
         infile: healpix format
                 LVC probability sky map
         percentage: float
                  probability percentage of the enclosed area  
         short_name: str
                 output file name
     
    Output: fits format
                 MOC map named "short_name"_"percentage" 
                 
                 Remark: for json format change the statement
                 "moc.write(short_name+'_MOC_'+str(percentage), format='fits' )" -->  
                "moc.write(short_name+'_MOC_'+str(percentage), format='json' )"        
    """
 
    import healpy as hp
    import numpy as np
    from math import log
    from astropy.table import Table
     
    #reading skymap
    hpx = hp.read_map( infile, verbose = False )
    npix = len( hpx )
    nside = hp.npix2nside( npix )
 
    sort = sorted( hpx, reverse = True )
    cumsum = np.cumsum( sort )
    index, value = min( enumerate( cumsum ), key = lambda x: abs( x[1] - percentage ) )

    # finding ipix indices confined in a given percentage 
    index_hpx = range( 0, len( hpx ) )
    hpx_index = np.c_[ hpx, index_hpx ]

    sort_2array = sorted( hpx_index, key = lambda x: x[0], reverse = True )
    value_contour = sort_2array[ 0:index ]

    j = 1 
    table_ipix_contour = [ ]

    for i in range ( 0, len( value_contour ) ):
        ipix_contour = int( value_contour[i][j] )
        table_ipix_contour.append( ipix_contour )
             
    # from index to polar coordinates
    theta, phi = hp.pix2ang( nside, table_ipix_contour )

    # converting these to right ascension and declination in degrees
    ra = np.rad2deg( phi )
    dec = np.rad2deg( 0.5 * np.pi - theta )

    # creating an astropy.table with RA[deg] and DEC[deg] ipix positions    
    contour_ipix = Table([ ra, dec ], names = ('RA[deg]', 'DEC[deg]'), 
                         meta = {'ipix': 'ipix table'})
        
    # setting MOC order    
    moc_order = int( log( nside, 2 ) )

    # creating a MOC map from the contour_ipix table
    moc = MOC.from_table( contour_ipix, 'RA[deg]', 'DEC[deg]', moc_order )

    # writing MOC file in fits
    moc.write( short_name + '_MOC_' + str( percentage ), format = 'fits' )

    # sending to Aladin plane
    aladin.send_file( short_name + '_MOC_' + str( percentage ) )
    aladin.rename ( plane = short_name  + '_MOC_' + str( percentage ) )

See the completed and updated results in the section Sky map visualization with Aladin

Questions, comments, suggestions, corrections, etc: email [email protected]

M. Branchesi, G. Greco and G. Stratta are supported by the Italian Ministry of Education, University and Research via grant FIRB 2012- RBFR12PM1F.