Sky Map I/O (ligo.skymap.io.fits)

Reading and writing HEALPix FITS files.

An example FITS header looks like this:

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

# HDU 1 in test.fits.gz
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                 4096 / length of dimension 1
NAXIS2  =                  192 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    1 / number of table fields
TTYPE1  = 'PROB    '
TFORM1  = '1024E   '
TUNIT1  = 'pix-1   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'RING    '           / Pixel ordering scheme, either RING or NESTED
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
EXTNAME = 'xtension'           / name of this binary table extension
NSIDE   =                  128 / Resolution parameter of HEALPIX
FIRSTPIX=                    0 / First pixel # (0 based)
LASTPIX =               196607 / Last pixel # (0 based)
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
OBJECT  = 'FOOBAR 12345'       / Unique identifier for this event
REFERENC= 'http://www.youtube.com/watch?v=0ccKPSVQcFk' / URL of this event
DATE-OBS= '2013-04-08T21:37:32.25' / UTC date of the observation
MJD-OBS =      56391.151064815 / modified Julian date of the observation
DATE    = '2013-04-08T21:50:32' / UTC date of file creation
CREATOR = 'fits.py '           / Program that created this file
RUNTIME =                 21.5 / Runtime in seconds of the CREATOR program
ligo.skymap.io.fits.read_sky_map(filename, nest=False, distances=False, moc=False, **kwargs)[source] [edit on github]

Read a LIGO/Virgo-type sky map and return a tuple of the HEALPix array and a dictionary of metadata from the header.

Parameters
filename: string

Path to the optionally gzip-compressed FITS file.

nest: bool, optional

If omitted or False, then detect the pixel ordering in the FITS file and rearrange if necessary to RING indexing before returning.

If True, then detect the pixel ordering and rearrange if necessary to NESTED indexing before returning.

If None, then preserve the ordering from the FITS file.

Regardless of the value of this option, the ordering used in the FITS file is indicated as the value of the ‘nest’ key in the metadata dictionary.

distances: bool, optional

If true, then read also read the additional HEALPix layers representing the conditional mean and standard deviation of distance as a function of sky location.

moc: bool, optional

If true, then preserve multi-order structure if present.

Examples

Test that we can read a legacy IDL-compatible file (https://bugs.ligo.org/redmine/issues/5168):

>>> import tempfile
>>> with tempfile.NamedTemporaryFile(suffix='.fits') as f:
...     nside = 512
...     npix = ah.nside_to_npix(nside)
...     ipix_nest = np.arange(npix)
...     hp.write_map(f.name, ipix_nest, nest=True, column_names=['PROB'])
...     m, meta = read_sky_map(f.name)
...     np.testing.assert_array_equal(m, hp.ring2nest(nside, ipix_nest))
ligo.skymap.io.fits.write_sky_map(filename, m, **kwargs)[source] [edit on github]

Write a gravitational-wave sky map to a file, populating the header with optional metadata.

Parameters
filename: str

Path to the optionally gzip-compressed FITS file.

mastropy.table.Table, numpy.array

If a Numpy record array or astorpy.table.Table instance, and has a column named ‘UNIQ’, then interpret the input as NUNIQ-style multi-order map [1]. Otherwise, interpret as a NESTED or RING ordered map.

**kwargs

Additional metadata to add to FITS header. If m is an astropy.table.Table instance, then the header is initialized from both m.meta and kwargs.

References

1

Górski, K.M., Wandelt, B.D., Hivon, E., Hansen, F.K., & Banday, A.J. 2017. The HEALPix Primer. The Unique Identifier scheme. http://healpix.sourceforge.net/html/intronode4.htm#SECTION00042000000000000000

Examples

Test header contents:

>>> order = 9
>>> nside = 2 ** order
>>> npix = ah.nside_to_npix(nside)
>>> prob = np.ones(npix, dtype=np.float) / npix
>>> import tempfile
>>> from ligo.skymap import version
>>> with tempfile.NamedTemporaryFile(suffix='.fits') as f:
...     write_sky_map(f.name, prob, nest=True,
...                   vcs_version='foo 1.0', vcs_revision='bar',
...                   build_date='2018-01-01T00:00:00')
...     for card in fits.getheader(f.name, 1).cards:
...         print(str(card).rstrip())
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                    8 / length of dimension 1
NAXIS2  =              3145728 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    1 / number of table fields
TTYPE1  = 'PROB    '
TFORM1  = 'D       '
TUNIT1  = 'pix-1   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NESTED  '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
NSIDE   =                  512 / Resolution parameter of HEALPIX
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
VCSVERS = 'foo 1.0 '           / Software version
VCSREV  = 'bar     '           / Software revision (Git)
DATE-BLD= '2018-01-01T00:00:00' / Software build date
>>> uniq = moc.nest2uniq(np.uint8(order), np.arange(npix))
>>> probdensity = prob / hp.nside2pixarea(nside)
>>> moc_data = np.rec.fromarrays(
...     [uniq, probdensity], names=['UNIQ', 'PROBDENSITY'])
>>> with tempfile.NamedTemporaryFile(suffix='.fits') as f:
...     write_sky_map(f.name, moc_data,
...                   vcs_version='foo 1.0', vcs_revision='bar',
...                   build_date='2018-01-01T00:00:00')
...     for card in fits.getheader(f.name, 1).cards:
...         print(str(card).rstrip())
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   16 / length of dimension 1
NAXIS2  =              3145728 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    2 / number of table fields
TTYPE1  = 'UNIQ    '
TFORM1  = 'K       '
TTYPE2  = 'PROBDENSITY'
TFORM2  = 'D       '
TUNIT2  = 'sr-1    '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NUNIQ   '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
MOCORDER=                    9 / MOC resolution (best order)
INDXSCHM= 'EXPLICIT'           / Indexing: IMPLICIT or EXPLICIT
VCSVERS = 'foo 1.0 '           / Software version
VCSREV  = 'bar     '           / Software revision (Git)
DATE-BLD= '2018-01-01T00:00:00' / Software build date