Posterior Samples to FITS Files (ligo-skymap-from-samples)

Example

This example shows how to extract the conditional distance posterior for a given sky location from the kernel density estimator.

First, run ligo-skymap-from-samples to create the KDE and the 3D sky map from the posterior sample chain:

$ ligo-skymap-from-samples --maxpts 1000 posterior_samples.dat

Then you can load the pickle in a Python interpreter and evaluate the KDE at any point:

>>> from astropy.coordinates import SkyCoord
>>> from matplotlib import pyplot as plt
>>> import numpy as np
>>> import pickle
>>> with open('skypost.obj', 'rb') as f:  
...     skypost = pickle.load(f)  
...
>>> coord = SkyCoord.from_name('NGC 4993')
>>> distance = np.arange(1, 100)
>>> coords = np.column_stack((np.tile(coord.ra.rad, len(distance)),
...                           np.tile(coord.dec.rad, len(distance)),
...                           distance))  
>>> post = skypost.posterior_spherical(coords)  
>>> plt.plot(distance, post)  
>>> plt.show()  

Generate a FITS sky map file from posterior samples using clustering and kernel density estimation.

The input file should be an HDF5 file with the following columns:

  • ra, rightascension, or right_ascension: J2000 right ascension in

    radians

  • dec or declination: J200 declination in radians

  • dist, distance, or luminosity_distance: luminosity distance in Mpc (optional)

The output consist of two files:

  • skypost.obj, a pickle representation of the kernel density estimator

  • skymap.fits.gz, a 3D localization in HEALPix/FITS format

usage: ligo-skymap-from-samples [-h] [--seed SEED] [--version]
                                [-l CRITICAL|ERROR|WARNING|INFO|DEBUG|NOTSET]
                                [--outdir OUTDIR]
                                [--fitsoutname SKYMAP.fits[.gz]]
                                [--loadpost SKYPOST.obj] [--maxpts MAXPTS]
                                [--trials TRIALS] [--enable-distance-map]
                                [--enable-multiresolution]
                                [--top-nside TOP_NSIDE] [-j [JOBS]]
                                [--instruments H1|L1|V1|... [H1|L1|V1|... ...]]
                                [--objid OBJID] [--path PATH]
                                [--tablename TABLENAME]
                                SAMPLES.hdf5

Positional Arguments

SAMPLES.hdf5

posterior samples file

Named Arguments

--version

show program’s version number and exit

-l, --loglevel

Default: INFO

--outdir, -o

output directory

Default: .

--fitsoutname

filename for the FITS file

Default: “skymap.fits”

--loadpost

filename for pickled posterior state

--maxpts

maximum number of posterior points to use; if omitted or greater than or equal to the number of posterior samples, then use all samples

--trials

number of trials at each clustering number

Default: 5

--enable-distance-map, --disable-distance-map

generate HEALPix map of distance estimates

Default: True

--enable-multiresolution, --disable-multiresolution

generate a multiresolution HEALPix map

Default: True

--top-nside

choose a start nside before HEALPix refinement steps (must be a valid nside)

Default: 16

-j, --jobs

Number of threads

Default: 1

--instruments

instruments to store in FITS header

--objid

event ID to store in FITS header

--path

The path of the dataset within the HDF5 file

--tablename

The name of the table to search for recursively within the HDF5 file. By default, search for posterior_samples

Default: “posterior_samples”

random number generator options

Options that affect the Numpy pseudo-random number genrator

--seed

Pseudo-random number generator seed [default: initialized from /dev/urandom or clock]