# Working with Multi-Order Sky Maps¶

For all events, LIGO/Virgo/KAGRA distributes both the standard HEALPix
format with the file extension `.fits.gz`

, as well as the new
multi-resolution HEALPix format, distinguished by the file extension
`.multiorder.fits`

. The multi-resolution HEALPix format is the primary and
preferred format, and the only format that is explicitly listed in the GCN
Notices and Circulars.

## What Problem Do Multi-Resolution Sky Maps Solve?¶

The multi-resolution format has been introduced as a forward-looking solution to deal with computational challenges related to highly accurate localizations. We are getting better at pinpointing gravitational-wave sources as more detectors come online and existing detectors become more sensitive. Unfortunately, as position accuracy improves, the size of the standard HEALPix sky maps will blow up. This started being a minor inconvenience in O2 with GW170817. It will get slowly worse as we approach design sensitivity. It’s already a major pain if you are studying future detector networks with simulations.

It is worth reviewing why LIGO/Virgo/KAGRA has adopted HEALPix rather than a more commonplace image format for sky maps in the first place. Gravitational-wave localizations are distinguished from many other kinds of astronomical image data by the following features:

The probability regions can subtend large angles.

They can wrap around the whole sky.

They can have multiple widely separated modes.

They can have irregular shapes or interference-like fringes.

As a consequence of these features, it is difficult to pick a good partial-sky projection (e.g. gnomonic, orthographic) in the general case. Traditional all-sky projections have wild variations in pixel size (e.g. plate carée) or shape (e.g. Mollweide, Aitoff) and as well as seams at the projection boundaries. HEALPix was already well-established for specialized uses in astronomy that cannot tolerate such projection artifacts (e.g. cosmic microwave background data sets and full-sky mosaics from optical surveys).

The natural sky resolution varies from one gravitational-wave event to another
depending on its SNR and the number of detectors. During early Advanced
LIGO and Virgo, HEALPix resolutions of `nside=512`

(\(6.9'\) per pixel)
to `nside=2048`

(\(1.7''\) per pixel) were adequate. However, as existing
gravitational-wave detectors improve in sensitivity and additional detectors
come online, finer resolutions will be required.

Fortunately, the increased resolution will come at little to no computational cost for actually producing localizations because most LIGO/Virgo/KAGRA parameter estimation analyses use a simple multi-resolution adaptive mesh refinement scheme that limits them to sampling the sky at only about 20k points.

When these multi-resolution meshes are flattened to a single HEALPix resolution, all but the finest nodes in the mesh become long sequences of repeated pixel values. High resolution also does not cost much in terms of disk space because gzip compression can store the long runs of repeated pixel values efficiently. The diagram below illustrates a multi-resolution structure that is fairly typical of gravitational-wave localizations from O1 and O2.

However, the resolution does come at a significant cost in the time it takes to decompress and read the FITS files (already up to tens of seconds for GW170817) and in terms of memory (up to several gigabytes). The time and memory will worsen as localization accuracy improves.

The multi-resolution format is immune to these issues because it is a direct representation of the adaptive mesh produced by the LIGO/Virgo/KAGRA localization algorithms.

## The UNIQ Indexing Scheme¶

Recall from before that three pieces of information are required to specify a HEALPix tile: nside to specify the resolution, ipix to identify a sky position at that resolution, and the indexing scheme.

HEALPix has a couple different indexing schemes. In the **RING** scheme,
indices advance west to east and then north to south. In the **NESTED** scheme,
indices encode the hierarchy of parent pixels in successively lower
resolutions. The image below illustrates these two indexing schemes.

There is a third HEALPix indexing scheme called **UNIQ**. The **UNIQ** indexing
scheme is special because it encodes both the resolution *and* the sky position
in a single integer. It assigns a single unique integer to every HEALPix tile
at every resolution. If ipix is the pixel index in the NESTED ordering, then
the unique pixel index uniq is:

The inverse is:

## FITS Format for Multi-Order Sky Maps¶

The FITS format for LIGO/Virgo/KAGRA multi-resolution sky maps uses the
**UNIQ** indexing scheme and is a superset of the FITS serialization for
Multi-Order Coverage (MOC) maps specified by IVOA [4] as part of
the Hierarchical Progressive Survey (HiPS) capability [3], notably
used by Aladin for storing and display all-sky image mosaics.

Let’s download an example multi-order FITS file with curl:

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

Let’s look at the FITS header:

```
$ fitsheader bayestar.multiorder.fits
# HDU 0 in bayestar.multiorder.fits:
SIMPLE = T / conforms to FITS standard
BITPIX = 8 / array data type
NAXIS = 0 / number of array dimensions
EXTEND = T
# HDU 1 in bayestar.multiorder.fits:
XTENSION= 'BINTABLE' / binary table extension
BITPIX = 8 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 40 / length of dimension 1
NAXIS2 = 19200 / length of dimension 2
PCOUNT = 0 / number of group parameters
GCOUNT = 1 / number of groups
TFIELDS = 5 / number of table fields
TTYPE1 = 'UNIQ '
TFORM1 = 'K '
TTYPE2 = 'PROBDENSITY'
TFORM2 = 'D '
TUNIT2 = 'sr-1 '
TTYPE3 = 'DISTMU '
TFORM3 = 'D '
TUNIT3 = 'Mpc '
TTYPE4 = 'DISTSIGMA'
TFORM4 = 'D '
TUNIT4 = 'Mpc '
TTYPE5 = 'DISTNORM'
TFORM5 = 'D '
TUNIT5 = 'Mpc-2 '
PIXTYPE = 'HEALPIX ' / HEALPIX pixelisation
ORDERING= 'NUNIQ ' / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C ' / Ecliptic, Galactic or Celestial (equatorial)
MOCORDER= 11 / MOC resolution (best order)
INDXSCHM= 'EXPLICIT' / Indexing: IMPLICIT or EXPLICIT
OBJECT = 'MS181101ab' / Unique identifier for this event
REFERENC= 'https://example.org/superevents/MS181101ab/view/' / URL of this event
INSTRUME= 'H1,L1,V1' / Instruments that triggered this event
DATE-OBS= '2018-11-01T22:22:46.654437' / UTC date of the observation
MJD-OBS = 58423.93248442635 / modified Julian date of the observation
DATE = '2018-11-01T22:34:49.000000' / UTC date of file creation
CREATOR = 'BAYESTAR' / Program that created this file
ORIGIN = 'LIGO/Virgo/KAGRA' / Organization responsible for this FITS file
RUNTIME = 3.24746292643249 / Runtime in seconds of the CREATOR program
DISTMEAN= 39.76999609489013 / Posterior mean distance (Mpc)
DISTSTD = 8.308435058808886 / Posterior standard deviation of distance (Mpc)
LOGBCI = 13.64819688928804 / Log Bayes factor: coherent vs. incoherent
LOGBSN = 261.0250944470225 / Log Bayes factor: signal vs. noise
VCSVERS = 'ligo.skymap 0.1.8' / Software version
VCSREV = 'becb07110491d799b753858845b5c24c82705404' / Software revision (Git)
DATE-BLD= '2019-07-25T22:36:58' / 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, min_inclination=0.0, max_inclination=1.5707963267948966, min_distance
HISTORY =None, max_distance=None, prior_distance_power=2, cosmology=False, mcmc=
HISTORY False, chain_dump=None, enable_snr_series=True, f_high_truncate=0.95)
HISTORY
HISTORY This was the command line that started the program:
HISTORY bayestar-localize-lvalert -N G298107 -o bayestar.multiorder.fits
```

This should look very similar to the FITS header for the standard HEALPix file from the previous section. The key differences are:

The

`ORDERING`

key has changed from`NESTED`

to`NUNIQ`

.The

`INDXSCHM`

key has changed from`IMPLICIT`

to`EXPLICIT`

.There is an extra column,

`UNIQ`

, that explicitly identifies each pixel in the**UNIQ**indexing scheme.The

`PROB`

column has been renamed to`PROBDENSITY`

, and the units have change from probability to probability per steradian.

## Reading Multi-Resolution Sky Maps¶

Now let’s go through some of the same common HEALPix operations from the
previous section, but using the multi-resolution format. Instead of Healpy, we
will use astropy-healpix because it has basic support for the **UNIQ**
indexing scheme.

First, we need the following imports:

```
>>> from astropy.table import QTable
>>> from astropy import units as u
>>> import astropy_healpix as ah
>>> import numpy as np
```

Next, let’s read the sky map. Instead of a special-purpose HEALPix method, we just read the FITS file into an Astropy table using Astropy’s unified file read/write interface:

```
>>> skymap = QTable.read('bayestar.multiorder.fits')
```

## Most Probable Sky Location¶

Next, let’s find the highest probability density sky position. This is a three-step process.

Find the

**UNIQ**pixel index of the highest probability density tile:`>>> i = np.argmax(skymap['PROBDENSITY']) >>> uniq = skymap[i]['UNIQ']`

What is the probability density per square degree in that tile?

`>>> skymap[i]['PROBDENSITY'].to_value(u.deg**-2) 0.07825164701914111`

Unpack the

**UNIQ**pixel index into the resolution,`nside`

, and the**NESTED**pixel index,`ipix`

, using the method`astropy_healpix.uniq_to_level_ipix()`

. (Note that this method returns`level`

, which is the logarithm base 2 of`nside`

, so we must also convert from`level`

to`nside`

using`astropy_healpix.level_to_nside()`

.)`>>> level, ipix = ah.uniq_to_level_ipix(uniq) >>> nside = ah.level_to_nside(level)`

Convert from

`nside`

and`ipix`

to right ascension and declination using`astropy_healpix.healpix_to_lonlat()`

(which is equivalent to`hp.pix2ang`

):`>>> ra, dec = ah.healpix_to_lonlat(ipix, nside, order='nested') >>> ra.deg 194.30419921874997 >>> dec.deg -17.856895095545468`

## Probability Density at a Known Position¶

Now let’s look up the probability density at a known sky position. In this case, let’s use the position of NGC 4993:

```
>>> ra = 197.4133 * u.deg
>>> dec = -23.3996 * u.deg
```

Brute Force Linear Search

The following brute force method of looking up a pixel by sky position has a complexity of \(O(N)\), where \(N\) is the number of multi-resolution pixels.

Unpack the

**UNIQ**pixel indices into their resolution and their**NESTED**pixel index.`>>> level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ']) >>> nside = ah.level_to_nside(level)`

Determine the

**NESTED**pixel index of the target sky position at the resolution of each multi-resolution tile.`>>> match_ipix = ah.lonlat_to_healpix(ra, dec, nside, order='nested')`

Find the multi-resolution tile whose

**NESTED**pixel index equals the target pixel index.`>>> i = np.flatnonzero(ipix == match_ipix)[0] >>> i 13484`

That pixel contains the target sky position. The probability density per square degree at that position is

`>>> skymap[i]['PROBDENSITY'].to_value(u.deg**-2) 0.034679190989078075`

Fast Binary Search

The following binary search method of looking up a pixel by sky position exploits the algebraic properties of HEALPix. It has a complexity of \(O(\log N)\) where \(N\) is the number of multi-resolution pixels. It assumes that every sky position is mapped on to exactly one multi-resolution tile, which is true for LIGO/Virgo/KAGRA multi-resolution sky maps.

First, find the

**NESTED**pixel index of every multi-resolution tile, at an arbitrarily high resolution. (`nside = 2**29`

works nicely because it is the highest possible HEALPix resolution that can be represented in a 64-bit signed integer.)`>>> max_level = 29 >>> max_nside = ah.level_to_nside(max_level) >>> level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ']) >>> index = ipix * (2**(max_level - level))**2`

Sort the pixels by this value.

`>>> sorter = np.argsort(index)`

Determine the

**NESTED**pixel index of the target sky location at that resolution.`>>> match_ipix = ah.lonlat_to_healpix(ra, dec, max_nside, order='nested')`

Do a binary search for that value.

`>>> i = sorter[np.searchsorted(index, match_ipix, side='right', sorter=sorter) - 1] >>> i 13484`

That pixel contains the target sky position. The probability density per square degree at that position is

`>>> skymap[i]['PROBDENSITY'].to_value(u.deg**-2) 0.034679190989078075`

## Find the 90% Probability Region¶

Follow these steps to find the 90% probability region.

Sort the pixels of the sky map by descending probability density.

`>>> skymap.sort('PROBDENSITY', reverse=True)`

Find the area of each pixel.

`>>> level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ']) >>> pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level))`

Calculate the probability within each pixel: the pixel area times the probability density.

`>>> prob = pixel_area * skymap['PROBDENSITY']`

Calculate the cumulative sum of the probability.

`>>> cumprob = np.cumsum(prob)`

Find the pixel for which the probability sums to 0.9 (90%).

`>>> i = cumprob.searchsorted(0.9)`

The area of the 90% credible region is simply the sum of the areas of the pixels up to that one.

`>>> area_90 = pixel_area[:i].sum() >>> area_90.to_value(u.deg**2) 30.975181093574633`

## Save the 90% Probability Region Footprint as a MOC¶

Having found the 90% probability region, it’s easy to now write it out its footprint as a MOC suitable for visualization in Aladin.

Keep only the pixels that are within the 90% credible region.

`>>> skymap = skymap[:i]`

Sort the pixels by their UNIQ pixel index.

`>>> skymap.sort('UNIQ')`

Delete all columns except for the UNIQ column.

`>>> skymap = skymap['UNIQ',]`

Write to a FITS file.

`>>> skymap.write('90percent.moc.fits')`