Merge pull request #18 from fermiPy/dmcat-dev

Dmcat dev

dmsky/skymap.py

@@ -7,10 +7,11 @@
 import healpy as hp
 import numpy as np
 
-from dmsky.utils.healpix import pix2ang, query_disc
+from dmsky.utils.healpix import pix2ang, query_disc, nside_to_order
 from dmsky.utils.coords import gal2cel
 from dmsky.roster import Roster
 
+from astropy.io import fits
 
 class Skymap(object):
     """Object to fill a HEALPix sky-map with J values from a `Roster`
@@ -36,14 +37,28 @@ def __init__(self, *targets, **kwargs):
             Coordinate system 'cel' or 'gal' ['gal']
 
         """
+        self.filename = kwargs.get('filename', None)
+        if self.filename is None:
+            self._construct(*targets, **kwargs)
+        else:
+            self._read_fits(filename)
+
+
+    def _construct(self, *targets, **kwargs):        
         self.nside = kwargs.get('nside', 512)
         self.coord = kwargs.get('coord', 'gal')
+        self.ann = kwargs.get('ann', True)
+        if self.ann:
+            self.units = 'GeV^2 cm^5 s^-2'
+        else:
+            self.units = 'GeV cm^3'
         self.values = np.zeros(hp.nside2npix(self.nside))
         self.pix = np.arange(hp.nside2npix(self.nside))
         self.lon, self.lat = pix2ang(self.nside, self.pix)
         self.roster = Roster(*targets)
         self._fill()
 
+
     def _fill(self):
         """Fills map values """
         for target in self.roster.values():
@@ -55,4 +70,50 @@ def _fill(self):
                 idx = query_disc(self.nside, target.ra, target.dec, target.psi_max)
                 lon, lat = self.lon[idx], self.lat[idx]
 
-            self.values[idx] += target.jvalue(lon, lat)
+            if self.ann:
+                self.values[idx] += target.jvalue(lon, lat)
+            else:
+                self.values[idx] += target.dvalue(lon, lat)
+
+
+    def write_fits(self, filename, overwrite=True):
+        """Write this Skymap to a FITS file"""
+        col = fits.Column(name='Bin0', format='1D', array=self.values)
+        header = fits.Header()
+        header['COORDSYS'] = self.coord
+        header['COORDTYPE'] = self.coord
+        header['ORDERING'] = 'RING'
+        header['INDXSCHM'] = 'IMPLICIT'
+        header['NSIDE'] = self.nside
+        header['ORDER'] = nside_to_order(self.nside)
+        header['FIRSTPIX'] = 0
+        header['LASTPIX'] = len(self.values) - 1
+        header['HPX_CONV'] = 'GALPROP'
+        header['BANDSHDU'] = ''
+
+        primhdu = fits.PrimaryHDU()
+        skyhdu = fits.BinTableHDU.from_columns([col], header=header, name='SKYMAP')
+
+        hdulist = fits.HDUList([primhdu, skyhdu])
+        hdulist.writeto(filename, overwrite=overwrite)
+
+    def _assert_header_key(self, header, key, value):
+        """Internal method to verify that a FITS header key matches the expected value"""
+        check = header.get(key)
+        if check != value:
+            raise ValueError("Header card %s=%s it should be %s" % (key, check, value))
+
+    def _read_fits(self, filename):
+        """Internal method to read this from a FITS file"""
+        f = fits.open(filename)
+        header = f[1].header
+        self._assert_header_key(header, 'HPX_CONV', 'GALPROP')
+        self._assert_header_key(header, 'ORDERING', 'RING')
+        self._assert_header_key(header, 'INDXSCHM', 'IMPLICIT')
+
+        self.nside = header['NSIDE']
+        self.coord = header['COORDSYS']
+        self.values = np.zeros(hp.nside2npix(self.nside))
+        self.pix = np.arange(hp.nside2npix(self.nside))
+        self.lon, self.lat = pix2ang(self.nside, self.pix)
+        self.values = f[2].data.field('Bin0')

dmsky/targets.py

@@ -183,7 +183,10 @@ def _d_integ(self):
     def _d_sigma(self):
         """Compute the uncertaintiy on the integrated D-factor
         """
-        dd = self.d_derivs
+        try:
+            dd = self.d_derivs
+        except ValueError:
+            return None
         den = self.density
         dv = np.matrix(np.zeros((len(den.deriv_params))))
         for i, pname in enumerate(den.deriv_params):
@@ -625,10 +628,14 @@ def write_jmap_wcs(self, filename, npix=150, clobber=False,
         self.setp('j_map_file', value=filename)
         return hdu.writeto(filename, clobber=clobber, **file_kwargs)
 
-    #def write_jmap_hpx(self, filename):
-    #    """Write the J-factor to a template map.
-    #    """
-    #    raise RuntimeError('write_jmap_hpx not implemented')
+    def write_jmap_hpx(self, filename):
+        """Write the J-factor to a template map.
+        """
+        from dmsky.skymap import Skymap
+        themap = Skymap([self], ann=True)
+        themap.write_fits(filename)
+        return themap
+
 
     write_jmap = write_jmap_wcs
 
@@ -663,10 +670,13 @@ def write_dmap_wcs(self, filename, npix=150, clobber=False,
         self.setp('d_map_file', value=filename)
         return hdu.writeto(filename, clobber=clobber, **file_kwargs)
 
-    #def write_dmap_hpx(self, filename):
-    #    """Write the D-factor to a template map.
-    #    """
-    #    raise RuntimeError('write_dmap_hpx not implemented')
+    def write_dmap_hpx(self, filename):
+        """Write the D-factor to a template map.
+        """
+        from dmsky.skymap import Skymap
+        themap = Skymap([self], ann=False)
+        themap.write_fits(filename)
+
 
     write_dmap = write_dmap_wcs
 

dmsky/utils/healpix.py

@@ -34,6 +34,10 @@ def subpixel(superpix, nside_superpix, nside_subpix):
 
 ############################################################
 
+NSIDE_TO_ORDER = {1:0, 2:1, 4:2, 8:3, 16:4, 32:5, 64:6, 128:7, 256:8, 512:9, 1024:10, 2048:11, 4096:12, 8192:13}
+
+def nside_to_order(nside):
+    return NSIDE_TO_ORDER.get(nside, -1)
 
 def pix2ang(nside, pix):
     """