Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org on Tuesday 22nd September 2020 starting at approximately 9am MST.It is expected to take around 15 minutes and there will be a short period of downtime towards the end of the maintenance window. Please address any comments, questions, or concerns to computing-help@igwn.org.

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)