Source code for gammapy.maps.hpxmap

# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import abc
import json
import numpy as np
from astropy.io import fits
from .base import Map
from .hpx import HpxGeom, HpxConv
from .utils import find_bintable_hdu, find_bands_hdu

__all__ = [
    'HpxMap',
]


[docs]class HpxMap(Map): """Base class for HEALPIX map classes. Parameters ---------- geom : `~gammapy.maps.HpxGeom` HEALPix geometry object. data : `~numpy.ndarray` Data array. meta : `~collections.OrderedDict` Dictionary to store meta data. """ def __init__(self, geom, data, meta=None): super(HpxMap, self).__init__(geom, data, meta) self._wcs2d = None self._hpx2wcs = None
[docs] @classmethod def create(cls, nside=None, binsz=None, nest=True, map_type='hpx', coordsys='CEL', data=None, skydir=None, width=None, dtype='float32', region=None, axes=None, conv='gadf', meta=None): """Factory method to create an empty HEALPix map. Parameters ---------- nside : int or `~numpy.ndarray` HEALPix NSIDE parameter. This parameter sets the size of the spatial pixels in the map. binsz : float or `~numpy.ndarray` Approximate pixel size in degrees. An NSIDE will be chosen that correponds to a pixel size closest to this value. This option is superseded by nside. nest : bool True for HEALPix "NESTED" indexing scheme, False for "RING" scheme. coordsys : {'CEL', 'GAL'}, optional Coordinate system, either Galactic ('GAL') or Equatorial ('CEL'). skydir : tuple or `~astropy.coordinates.SkyCoord` Sky position of map center. Can be either a SkyCoord object or a tuple of longitude and latitude in deg in the coordinate system of the map. map_type : {'hpx', 'hpx-sparse'} Map type. Selects the class that will be used to instantiate the map. width : float Diameter of the map in degrees. If None then an all-sky geometry will be created. axes : list List of `~MapAxis` objects for each non-spatial dimension. conv : {'fgst-ccube','fgst-template','gadf'}, optional Default FITS format convention that will be used when writing this map to a file. Default is 'gadf'. meta : `~collections.OrderedDict` Dictionary to store meta data. Returns ------- map : `~HpxMap` A HPX map object. """ from .hpxnd import HpxNDMap from .hpxsparse import HpxSparseMap hpx = HpxGeom.create(nside=nside, binsz=binsz, nest=nest, coordsys=coordsys, region=region, conv=conv, axes=axes, skydir=skydir, width=width) if cls.__name__ == 'HpxNDMap': return HpxNDMap(hpx, dtype=dtype, meta=meta) elif cls.__name__ == 'HpxSparseMap': return HpxSparseMap(hpx, dtype=dtype, meta=meta) elif map_type == 'hpx': return HpxNDMap(hpx, dtype=dtype, meta=meta) elif map_type == 'hpx-sparse': return HpxSparseMap(hpx, dtype=dtype, meta=meta) else: raise ValueError('Unrecognized map type: {}'.format(map_type))
[docs] @classmethod def from_hdulist(cls, hdu_list, hdu=None, hdu_bands=None): """Make a HpxMap object from a FITS HDUList. Parameters ---------- hdu_list : `~astropy.io.fits.HDUList` HDU list containing HDUs for map data and bands. hdu : str Name or index of the HDU with the map data. If None then the method will try to load map data from the first BinTableHDU in the file. hdu_bands : str Name or index of the HDU with the BANDS table. Returns ------- hpx_map : `HpxMap` Map object """ if hdu is None: hdu_out = find_bintable_hdu(hdu_list) else: hdu_out = hdu_list[hdu] if hdu_bands is None: hdu_bands = find_bands_hdu(hdu_list, hdu_out) hdu_bands_out = None if hdu_bands is not None: hdu_bands_out = hdu_list[hdu_bands] return cls.from_hdu(hdu_out, hdu_bands_out)
[docs] def to_hdulist(self, hdu='SKYMAP', hdu_bands=None, sparse=False, conv=None): """Convert to `~astropy.io.fits.HDUList`. Parameters ---------- hdu : str The HDU extension name. hdu_bands : str The HDU extension name for BANDS table. sparse : bool Set INDXSCHM to SPARSE and sparsify the map by only writing pixels with non-zero amplitude. conv : {'fgst-ccube','fgst-template','gadf',None}, optional FITS format convention. If None this will be set to the default convention of the map. Returns ------- hdu_list : `~astropy.io.fits.HDUList` """ if self.geom.axes: hdu_bands_out = self.geom.make_bands_hdu(hdu=hdu_bands, hdu_skymap=hdu, conv=conv) hdu_bands = hdu_bands_out.name else: hdu_bands_out = None hdu_bands = None hdu_out = self.make_hdu(hdu=hdu, hdu_bands=hdu_bands, sparse=sparse, conv=conv) hdu_out.header['META'] = json.dumps(self.meta) hdu_list = [fits.PrimaryHDU(), hdu_out] if self.geom.axes: hdu_list += [hdu_bands_out] return fits.HDUList(hdu_list)
[docs] @abc.abstractmethod def to_wcs(self, sum_bands=False, normalize=True, proj='AIT', oversample=2, width_pix=None, hpx2wcs=None): """Make a WCS object and convert HEALPIX data into WCS projection. Parameters ---------- sum_bands : bool Sum over non-spatial axes before reprojecting. If False then the WCS map will have the same dimensionality as the HEALPix one. normalize : bool True -> preserve integral by splitting HEALPIX values between bins proj : str WCS-projection oversample : float Oversampling factor for WCS map. This will be the approximate ratio of the width of a HPX pixel to a WCS pixel. If this parameter is None then the width will be set from ``width_pix``. width_pix : int Width of the WCS geometry in pixels. The pixel size will be set to the number of pixels satisfying ``oversample`` or ``width_pix`` whichever is smaller. If this parameter is None then the width will be set from ``oversample``. hpx2wcs : `~HpxToWcsMapping` Set the HPX to WCS mapping object that will be used to generate the WCS map. If none then a new mapping will be generated based on ``proj`` and ``oversample`` arguments. Returns ------- map_out : `~gammapy.maps.WcsMap` WCS map object. """ pass
[docs] @abc.abstractmethod def to_swapped(self): """Return a new map with the opposite scheme (ring or nested). Returns ------- map : `~HpxMap` Map object. """ pass
[docs] @abc.abstractmethod def to_ud_graded(self, nside, preserve_counts=False): """Upgrade or downgrade the resolution of the map to the chosen nside. Parameters ---------- nside : int NSIDE parameter of the new map. preserve_counts : bool Choose whether to preserve counts (total amplitude) or intensity (amplitude per unit solid angle). Returns ------- map : `~HpxMap` Map object. """ pass
[docs] def make_hdu(self, hdu=None, hdu_bands=None, sparse=False, conv=None): """Make a FITS HDU with input data. Parameters ---------- hdu : str The HDU extension name. hdu_bands : str The HDU extension name for BANDS table. sparse : bool Set INDXSCHM to SPARSE and sparsify the map by only writing pixels with non-zero amplitude. conv : {'fgst-ccube','fgst-template','gadf',None}, optional FITS format convention. If None this will be set to the default convention of the map. Returns ------- hdu_out : `~astropy.io.fits.BinTableHDU` or `~astropy.io.fits.ImageHDU` Output HDU containing map data. """ convname = self.geom.conv if conv is None else conv conv = HpxConv.create(convname) hduname = conv.hduname if hdu is None else hdu hduname_bands = conv.bands_hdu if hdu_bands is None else hdu_bands header = self.geom.make_header(conv=conv) if self.geom.axes: header['BANDSHDU'] = hduname_bands if sparse: header['INDXSCHM'] = 'SPARSE' cols = [] if header['INDXSCHM'] == 'EXPLICIT': array = self.geom._ipix cols.append(fits.Column('PIX', 'J', array=array)) elif header['INDXSCHM'] == 'LOCAL': array = np.arange(self.data.shape[-1]) cols.append(fits.Column('PIX', 'J', array=array)) cols += self._make_cols(header, conv) return fits.BinTableHDU.from_columns(cols, header=header, name=hduname)