-
Notifications
You must be signed in to change notification settings - Fork 1
/
artquicklook_healpix.py
57 lines (48 loc) · 1.68 KB
/
artquicklook_healpix.py
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
import numpy as np
import matplotlib.pyplot as plt
import os
import arttools
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
import healpy as hp
import os, shutil, argparse
parser = argparse.ArgumentParser()
parser.add_argument("--list", help="input list of L1b files")
parser.add_argument("--mapname", help="name of map to create, default 'artxc_survmap.fits.gz'", default='artxc_survmap.fits.gz')
args = parser.parse_args()
inputlist = args.list
if inputlist==None or not inputlist:
print ('>>ERROR>> Please, provide list of input L1b files')
exit(1)
else:
print ('>>>>>>>>> input list:'+inputlist)
mapname = args.mapname
if mapname==None or not mapname:
print ('>>>>>>>>> will create artxc_survmap.fits.gz')
else:
print ('>>>>>>>>> will create '+mapname)
inp = open(inputlist, 'r')
fnames = inp.readlines()
inp.close()
NSIDE = 1024*2
NPIX = hp.nside2npix(NSIDE)
sky = np.zeros(NPIX)
for fname in fnames[:]:
print('adding file: ',fname.strip())
urdfile = fits.open(fname.strip())
urddata = urdfile[1].data
energy = urddata['ENERGY']
grades = urddata['GRADE']
flags = urddata['FLAG']
mask = np.bitwise_and(energy>=4, energy<=11.)
gmask = np.bitwise_and(np.bitwise_and(grades>=0, grades<=9.), flags==0)
mask = np.bitwise_and(mask, gmask)
ra, dec = urddata['RA'][mask], urddata['DEC'][mask]
coords = SkyCoord(ra*u.degree, dec*u.degree, frame='icrs', unit='deg')
l ,b = coords.galactic.l.degree, coords.galactic.b.degree
pix = hp.ang2pix(NSIDE, l, b, lonlat=True)
for p in pix:
sky[p]+=1
hp.write_map(mapname, sky, overwrite=True)