Working with Sky Maps

Let’s take a look at what is inside one of the LIGO/Virgo probability sky maps. They are FITS image files and can be manipulated and viewed with many commonplace FITS tools. However, they are a little unusual in two regards. First, since they are all-sky images, they are stored in the HEALPix projection, a format that is used for Planck all-sky CMB maps and by Aladin for hierarchical, progressively refined, all-sky survey images (HiPS). Second, the value stored at each pixel is the probability that the gravitational-wave source is within that pixel.

HEALPix projection

Let’s download an example FITS file with curl:

$ curl -O https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz

We can look at the metadata inside the FITS file by printing its header with tools like funhead from Funtools, imhead from WCSTools, or fitsheader from Astropy:

$ fitsheader bayestar.fits.gz
# HDU 0 in bayestar.fits.gz:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                    8 / array data type
NAXIS   =                    0 / number of array dimensions
EXTEND  =                    T

# HDU 1 in bayestar.fits.gz:
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   32 / length of dimension 1
NAXIS2  =               786432 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    4 / number of table fields
TTYPE1  = 'PROB    '
TFORM1  = 'D       '
TUNIT1  = 'pix-1   '
TTYPE2  = 'DISTMU  '
TFORM2  = 'D       '
TUNIT2  = 'Mpc     '
TTYPE3  = 'DISTSIGMA'
TFORM3  = 'D       '
TUNIT3  = 'Mpc     '
TTYPE4  = 'DISTNORM'
TFORM4  = 'D       '
TUNIT4  = 'Mpc-2   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NESTED  '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
NSIDE   =                  256 / Resolution parameter of HEALPIX
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
OBJECT  = 'M2052   '           / Unique identifier for this event
REFERENC= 'https://gracedb-playground.ligo.org/events/M2052' / URL of this event
INSTRUME= 'H1,L1   '           / Instruments that triggered this event
DATE-OBS= '2018-11-01T22:22:46.654438' / UTC date of the observation
MJD-OBS =    58423.93248442614 / modified Julian date of the observation
DATE    = '2018-11-01T22:34:37.000000' / UTC date of file creation
CREATOR = 'BAYESTAR'           / Program that created this file
ORIGIN  = 'LIGO/Virgo'         / Organization responsible for this FITS file
RUNTIME =                 11.0 / Runtime in seconds of the CREATOR program
DISTMEAN=    141.1453950128411 / Posterior mean distance (Mpc)
DISTSTD =    39.09548411497191 / Posterior standard deviation of distance (Mpc)
LOGBCI  =    7.793862946657789 / Log Bayes factor: coherent vs. incoherent
LOGBSN  =    47.28194676827084 / Log Bayes factor: signal vs. noise
VCSVERS = 'ligo.skymap 0.0.17' / Software version
VCSREV  = 'cb59e5fd04d41c5181ae9e41fe59de232877ddd2' / Software revision (Git)
DATE-BLD= '2018-10-24T20:50:55' / Software build date
HISTORY
HISTORY Generated by calling the following Python function:
HISTORY ligo.skymap.bayestar.localize(event=..., waveform='o2-uberbank', f_low=3
HISTORY 0.0, min_distance=None, max_distance=None, prior_distance_power=None, co
HISTORY smology=False, mcmc=False, chain_dump=None, enable_snr_series=True, f_hi
HISTORY gh_truncate=0.95)
HISTORY
HISTORY This was the command line that started the program:
HISTORY gwcelery worker -l info -n gwcelery-openmp-worker -Q openmp -c 1

There are several useful pieces of information here:

  • COORDSYS=C, telling you that the HEALPix projection is in the Celestial (equatorial, ICRS) frame, as all LIGO/Virgo probability sky maps will be.
  • OBJECT, the unique LIGO/Virgo identifier for the event.
  • REFERENC, a link to the candidate page in GraceDb.
  • INSTRUME, a list of gravitational-wave sites that triggered on the event: H1 for LIGO Hanford, L1 for LIGO Livingston, and V1 for Virgo.
  • DATE-OBS, the UTC time of the event. In the case of a compact binary coalescence candidate, this is the time that the signal from the merger passed through the geocenter.
  • MJD-OBS, same as DATE-OBS, but given as a modified Julian day.

You can view the sky map in many common FITS image viewers such as Aladin:

Aladin screenshot

or DS9 (although DS9 shows HEALPix sky maps in an unusual orientation; see Figure 4 of Calabretta & Roukema (2007) for more information.

DS9 screenshot

Now, let’s go through some examples of manipulating HEALPix sky maps programmatically. The HEALPix project provides official libraries for many languages, including C, C++, Fortran, IDL, and Java. However, since this is a Python tutorial, we are going to demonstrate how to manipulate HEALPix maps with the official Python library, Healpy.

Reading Sky Maps

First, if you have not already downloaded an example sky map, you can do so now by having Python call curl on the command line:

Next, we need to read in the file in Python with Healpy:

>>> hpx = hp.read_map('bayestar.fits.gz')
NSIDE = 256
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING

You can suppress printing informational messages while loading the file by passing the keyword argument verbose=False. You can read both the HEALPix image data and the FITS header by passing the h=True keyword argument:

>>> hpx, header = hp.read_map('bayestar.fits.gz', h=True, verbose=False)

Manipulating HEALPix Coordinates

The image data is a 1D array of values:

>>> hpx
array([6.22405744e-25, 1.46981290e-25, 1.94449365e-25, ...,
       2.33147793e-20, 6.78207416e-21, 3.07118068e-22])

Healpy has several useful plotting routines including hp.mollview for plotting a Mollweide-projection all-sky map:

>>> hp.mollview(hpx)

Each entry in the array represents the probability contained within a quadrilateral pixel whose position on the sky is uniquely specified by the index in the array and the array’s length. Because HEALPix pixels are equal area, we can find the number of pixels per square degree just from the length of the HEALPix array:

>>> npix = len(hpx)
>>> sky_area = 4 * 180**2 / np.pi
>>> sky_area / npix
0.052455852825697924

The function hp.pix2ang converts from pixel index to spherical polar coordinates; the function hp.ang2pix does the reverse.

Both hp.pix2ang and hp.ang2pix take, as their first argument, nside, the lateral resolution fo the HEALPix map. You can find nside from the length of the image array by calling hp.npix2nside:

>>> nside = hp.npix2nside(npix)
>>> nside
256

Let’s look up the right ascension and declination of pixel number 123. We’ll call hp.pix2ang to get the spherical polar coordinates \((\theta, \phi)\) in radians, and then use np.rad2deg to convert these to right ascension and declination in degrees.

>>> ipix = 123
>>> theta, phi = hp.pix2ang(nside, ipix)
>>> ra = np.rad2deg(phi)
>>> dec = np.rad2deg(0.5 * np.pi - theta)
>>> ra, dec
(129.375, 88.5380288373519)

Let’s find which pixel contains the point RA=194.95, Dec=27.98.

>>> ra = 194.95
>>> dec = 27.98
>>> theta = 0.5 * np.pi - np.deg2rad(dec)
>>> phi = np.deg2rad(ra)
>>> ipix = hp.ang2pix(nside, theta, phi)
>>> ipix
208938

Most Probable Sky Location

Let’s find the highest probability pixel. What is the probability inside it?

>>> ipix_max = np.argmax(hpx)
>>> hpx[ipix_max]
9.35702310989353e-05

Where is the highest probability pixel on the sky? Use hp.pix2ang.

>>> theta, phi = hp.pix2ang(nside, ipix_max)
>>> ra = np.rad2deg(phi)
>>> dec = np.rad2deg(0.5 * np.pi - theta)
>>> ra, dec
(90.87890625, -40.620185190672686)

Integrated Probability in a Circle

How do we find the probability that the source is contained within a circle on the sky? First we find the pixels that are contained within the circle using hp.query_disc. Note that this function takes as its arguments the Cartesian coordinates of the center of the circle, and its radius in radians. Then, we sum the values of the HEALPix image array contained at those pixels.

First, we define the RA, Dec, and radius of circle in degrees:

>>> ra = 213.22
>>> dec = -37.45
>>> radius = 3.1

Then we convert to spherical polar coordinates and radius of circle in radians:

>>> theta = 0.5 * np.pi - np.deg2rad(dec)
>>> phi = np.deg2rad(ra)
>>> radius = np.deg2rad(radius)

Then we calculate the Cartesian coordinates of the center of circle:

>>> xyz = hp.ang2vec(theta, phi)

We call hp.query_disc, which returns an array of the indices of the pixels that are inside the circle:

>>> ipix_disc = hp.query_disc(nside, xyz, radius)

Finally, we sum the probability in all of the matching pixels:

>>> hpx[ipix_disc].sum()
9.522375325439142e-06

Integrated Probability in a Polygon

Similarly, we can use the hp.query_polygon function to look up the indices of the pixels within a polygon (defined by the Cartesian coordinates of its vertices), and then compute the probability that the source is inside that polygon by summing the values of the pixels.

>>> xyz = [[-0.69601758, -0.41315628, -0.58724902],
...        [-0.68590811, -0.40679797, -0.60336181],
...        [-0.69106913, -0.39820114, -0.60320752],
...        [-0.7011786 , -0.40455945, -0.58709473]]
>>> ipix_poly = hp.query_polygon(nside, xyz)
>>> hpx[ipix_poly].sum()
3.935524328237466e-11

These are all of the HEALPix functions from Healpy that we will need for the remainder of the this tutorial.

Other useful Healpy functions include hp.ud_grade for upsampling or downsampling a sky map and hp.get_interp_val for performing bilinear interpolation between pixels. See the Healpy tutorial for other useful operations.