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
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)
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  as part of the Hierarchical Progressive Survey (HiPS) capability , 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:
ORDERINGkey has changed from
INDXSCHMkey has changed from
There is an extra column,
UNIQ, that explicitly identifies each pixel in the UNIQ indexing scheme.
PROBcolumn 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 Table >>> from astropy import units as u >>> import astropy_healpix as ah >>> import numpy as np
>>> skymap = Table.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'] * (np.pi / 180)**2 0.0782516470191411
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, ipix = ah.uniq_to_level_ipix(uniq) >>> nside = ah.level_to_nside(level)
>>> 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) >>> i 13484
That pixel contains the target sky position.
>>> skymap[i]['PROBDENSITY'] * (np.pi / 180)**2 0.03467919098907807
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**29works 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.
>>> skymap[i]['PROBDENSITY'] * (np.pi / 180)**2 0.03467919098907807