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
(nside=2048
(
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 fromNESTED
toNUNIQ
.The
INDXSCHM
key has changed fromIMPLICIT
toEXPLICIT
.There is an extra column,
UNIQ
, that explicitly identifies each pixel in the UNIQ indexing scheme.The
PROB
column has been renamed toPROBDENSITY
, 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 methodastropy_healpix.uniq_to_level_ipix()
. (Note that this method returnslevel
, which is the logarithm base 2 ofnside
, so we must also convert fromlevel
tonside
usingastropy_healpix.level_to_nside()
.)>>> level, ipix = ah.uniq_to_level_ipix(uniq) >>> nside = ah.level_to_nside(level)
Convert from
nside
andipix
to right ascension and declination usingastropy_healpix.healpix_to_lonlat()
(which is equivalent tohp.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
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
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')