Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org starting at around 10am CST on Tuesday 28 January 2020. It is expected to take around an hour and there will be a short period, around five minutes, of downtime towards the end of the maintenance window. In addition the GitLab-CI runners will be paused starting at around 9am CST in order to be updated.

lalinference_average_skymaps.py 2.76 KB
Newer Older
Adam Mercer's avatar
Adam Mercer committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
#
# Copyright (C) 2017  Leo Singer
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
"""
Average several sky maps according to the Burst group established by:
 * LIGO-T1700050-v7 (https://dcc.ligo.org/LIGO-T1700050-v7)
 * LIGO-T1700078-v2 (https://dcc.ligo.org/LIGO-T1700078-v2)

The sky maps are first upsampled to a common resolution and then averaged
pixel-by-pixel. Some FITS header values are copied from the first input
sky map.
"""
__author__ = "Leo Singer <leo.singer@ligo.org>"


# Command line interface
import argparse
from lalinference.bayestar import command
parser = command.ArgumentParser()
parser.add_argument(
    'input', metavar='INPUT.fits[.gz]', type=argparse.FileType('rb'),
    nargs='+', help='Input FITS files')
parser.add_argument(
    '-o', '--output', metavar='OUTPUT.fits[.gz]',
    help='Ouptut FITS file [default: derived from input filenames]')
opts = parser.parse_args()


# Late imports
import os.path
import healpy as hp
import numpy as np
from lalinference.io import fits


# Read all input sky maps
probs, metas = zip(*(fits.read_sky_map(file.name, nest=True)
                     for file in opts.input))

# Determine output HEALPix resolution
npix = max(len(prob) for prob in probs)
nside = hp.npix2nside(npix)

# Average sky maps
prob = np.mean([hp.ud_grade(prob, nside,
                            order_in='NESTED', order_out='NESTED',
                            power=-2) for prob in probs], axis=0)

# Describe the processing of this file
history = ['Arithmetic mean of the following {} sky maps:'.format(len(probs))]
history += ['    ' + file.name for file in opts.input]

# Create metadata dictionary for FITS header
meta = dict(creator=parser.prog, nest=True, history=history)

# Copy some header values from the first sky map
copy_keys = ('objid', 'url', 'instruments', 'gps_time')
meta.update({key: metas[0][key] for key in copy_keys if key in metas[0]})

# Determine output filename
if opts.output is not None:
    output = opts.output
else:
    output = '_plus_'.join(os.path.split(file.name)[1].partition('.fits')[0]
                      for file in opts.input) + '.fits.gz'

# Write output filename to stdout
print(output)

# Write output file
fits.write_sky_map(output, prob, **meta)