gammapy-coordinate-images

#!/usr/bin/env python
# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Make maps that can be used to create profiles.

The following images can be created:
* LON -- Longitude coordinate
* LAT -- Latitude coordinate
* DIST -- Distance to mask
* SOLID_ANGLE -- Solid angle
"""

# Parse command line arguments

from gammapy.utils.scripts import argparse, GammapyFormatter
parser = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=GammapyFormatter)
parser.add_argument('infile', type=str,
                    help='Input FITS file name')
parser.add_argument('outfile', type=str,
                    help='Output FITS file name')
parser.add_argument('--make_coordinate_maps', action='store_true',
                    help='Create coordinate maps')
parser.add_argument('--make_distance_map', action='store_true',
                    help='Create distance to mask map')
parser.add_argument('--clobber', action='store_true',
                    help='Clobber output files?')

args = parser.parse_args()
args = vars(args)

# Execute script

import logging
logging.basicConfig(level=logging.DEBUG, format='%(levelname)s - %(message)s')
from astropy.io import fits
from gammapy.utils.fits import get_hdu

logging.info('Reading {0}'.format(args['infile']))
hdu = get_hdu(args['infile'])

out_hdus = fits.HDUList()

if args['make_coordinate_maps']:
    from gammapy.utils.image import coordinates
    logging.info('Computing LON and LAT maps')
    lon, lat = coordinates(hdu)
    out_hdus.append(fits.ImageHDU(lon, hdu.header, 'LON'))
    out_hdus.append(fits.ImageHDU(lat, hdu.header, 'LAT'))

if args['make_distance_map']:
    from gammapy.utils.image import exclusion_distance
    logging.info('Computing DIST map')
    dist = exclusion_distance(hdu.data)
    out_hdus.append(fits.ImageHDU(dist, hdu.header, 'DIST'))

logging.info('Writing {0}'.format(args['outfile']))
out_hdus.writeto(args['outfile'], clobber=True)