From 50ec027249485c6de2f9d60dec717e8f90419e33 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Tue, 9 Jul 2019 15:18:59 -0700 Subject: [PATCH 1/4] Added nside_to_order function --- dmsky/utils/healpix.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dmsky/utils/healpix.py b/dmsky/utils/healpix.py index b7bdc20..634b555 100644 --- a/dmsky/utils/healpix.py +++ b/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): """ From 2aee28ff8e5d3b03714bd991ec028bfeb6886a55 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Tue, 9 Jul 2019 15:19:31 -0700 Subject: [PATCH 2/4] Added write_jmap_hpx implemenation in targets.py --- dmsky/targets.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/dmsky/targets.py b/dmsky/targets.py index d710f86..354ecce 100644 --- a/dmsky/targets.py +++ b/dmsky/targets.py @@ -625,10 +625,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 +667,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 From 176349eb4f29259ac5d618df73caab1b84ef3b39 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Tue, 9 Jul 2019 15:19:56 -0700 Subject: [PATCH 3/4] Added functionality to Skymap --- dmsky/skymap.py | 65 +++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 63 insertions(+), 2 deletions(-) diff --git a/dmsky/skymap.py b/dmsky/skymap.py index 065efad..ba20967 100644 --- a/dmsky/skymap.py +++ b/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') From 87a67e074fbb5f620c217857f01b68e742155f21 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Thu, 3 Sep 2020 10:53:10 -0700 Subject: [PATCH 4/4] catch error in making d-factors --- dmsky/targets.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/dmsky/targets.py b/dmsky/targets.py index 354ecce..e5244c5 100644 --- a/dmsky/targets.py +++ b/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):