|
1
|
|
|
from hawc_hal.maptree import map_tree_factory |
|
2
|
|
|
from astropy.io import fits |
|
3
|
|
|
|
|
4
|
|
|
from IPython import embed |
|
5
|
|
|
|
|
6
|
|
|
import healpy as hp |
|
7
|
|
|
import numpy as np |
|
8
|
|
|
|
|
9
|
|
|
import argparse |
|
10
|
|
|
import os |
|
11
|
|
|
from datetime import datetime |
|
12
|
|
|
|
|
13
|
|
|
|
|
14
|
|
|
parser = argparse.ArgumentParser(description='Process some integers.') |
|
15
|
|
|
|
|
16
|
|
|
parser.add_argument('--input', metavar='input', type=str, |
|
17
|
|
|
required=True, |
|
18
|
|
|
help='pick an hdf5 file as input') |
|
19
|
|
|
|
|
20
|
|
|
parser.add_argument('--output', metavar='output', type=str, |
|
21
|
|
|
default="map", |
|
22
|
|
|
help='pick a file prefix (e.g. given "test" -> test_binX.fits.gz)') |
|
23
|
|
|
|
|
24
|
|
|
parser.add_argument('--overwrite', |
|
25
|
|
|
action="store_true", |
|
26
|
|
|
help='Overwrite existing files') |
|
27
|
|
|
|
|
28
|
|
|
|
|
29
|
|
|
args = parser.parse_args() |
|
30
|
|
|
|
|
31
|
|
|
filename=os.path.expandvars(args.input) |
|
32
|
|
|
outfile=os.path.expandvars(args.output) |
|
33
|
|
|
clobber=args.overwrite |
|
34
|
|
|
|
|
35
|
|
|
if not os.path.isfile(filename): |
|
36
|
|
|
print("You must enter a valid file!") |
|
37
|
|
|
exit(1) |
|
38
|
|
|
|
|
39
|
|
|
|
|
40
|
|
|
# Export the entire map tree (full sky) |
|
41
|
|
|
# this throws a warning for partial map trees but its OK |
|
42
|
|
|
maptree = map_tree_factory(filename, None) |
|
43
|
|
|
|
|
44
|
|
|
|
|
45
|
|
|
now=datetime.now() |
|
46
|
|
|
startMJD=56987.9286332 |
|
47
|
|
|
|
|
48
|
|
|
|
|
49
|
|
|
#FIRST HEADER |
|
50
|
|
|
''' |
|
51
|
|
|
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy |
|
52
|
|
|
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H |
|
53
|
|
|
DATE = '2018-12-01T02:31:14' / file creation date (YYYY-MM-DDThh:mm:ss UT) |
|
54
|
|
|
STARTMJD= 56987.9286332451 / MJD of first event |
|
55
|
|
|
STOPMJD = 58107.2396848326 / MJD of last event |
|
56
|
|
|
NEVENTS = -1. / Number of events in map |
|
57
|
|
|
TOTDUR = 24412.9020670185 / Total integration time [hours] |
|
58
|
|
|
DURATION= 1.9943578604616 / Avg integration time [hours] |
|
59
|
|
|
MAPTYPE = 'duration' / e.g. Skymap, Moonmap |
|
60
|
|
|
MAXDUR = -1. / Max integration time [hours] |
|
61
|
|
|
MINDUR = -1. / Min integration time [hours] |
|
62
|
|
|
EPOCH = 'unknown ' / e.g. J2000, current, J2016, B1950, etc. |
|
63
|
|
|
HIERARCH MAPFILETYPE = 'duration' / e.g. standard, duration, integration |
|
64
|
|
|
''' |
|
65
|
|
|
|
|
66
|
|
|
FITS_COMMENT="FITS (Flexible Image Transport System) format is defined in 'Astronomy and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H" |
|
67
|
|
|
|
|
68
|
|
|
primary_keys=['COMMENT', 'COMMENT', 'DATE', 'STARTMJD', 'STOPMJD', |
|
69
|
|
|
'NEVENTS', 'TOTDUR', 'DURATION', 'MAPTYPE', 'MAXDUR', |
|
70
|
|
|
'MINDUR', 'EPOCH', 'MAPFILETYPE'] |
|
71
|
|
|
|
|
72
|
|
|
primary_values=[FITS_COMMENT , FITS_COMMENT, "{0}".format(now), 56987.9286332451, 58107.2396848326, |
|
73
|
|
|
-1.0, 24412.9020670185, 1.9943578604616, 'duration', -1.0, -1.0, 'unknown', 'duration'] |
|
74
|
|
|
|
|
75
|
|
|
primary_comments=["file does conform to FITS standard", |
|
76
|
|
|
"number of bits per data pixel", |
|
77
|
|
|
"number of data axes", |
|
78
|
|
|
"FITS dataset may contain extension", |
|
79
|
|
|
"MJD of first event", "MJD of last event", |
|
80
|
|
|
"Number of events in map", |
|
81
|
|
|
"Total integration time [hours]", |
|
82
|
|
|
"Avg integration time [hours]", |
|
83
|
|
|
"e.g. Skymap, Moonmap", |
|
84
|
|
|
"Max integration time [hours]", |
|
85
|
|
|
"Min integration time [hours]", |
|
86
|
|
|
"e.g. J2000, current, J2016, B1950, etc.", |
|
87
|
|
|
"e.g. standard, duration, integration"] |
|
88
|
|
|
|
|
89
|
|
|
|
|
90
|
|
|
|
|
91
|
|
|
labels=['data map', 'background map', 'exposure map'] |
|
92
|
|
|
label_format=[ np.float64 for i in range(len(labels)) ] |
|
93
|
|
|
label_units=[ 'unknown' for i in range(len(labels)) ] |
|
94
|
|
|
|
|
95
|
|
|
|
|
96
|
|
|
for i, analysis_bin in enumerate(maptree.analysis_bins_labels): |
|
97
|
|
|
map_bin = maptree[analysis_bin] |
|
98
|
|
|
#properties |
|
99
|
|
|
nside = map_bin.nside |
|
100
|
|
|
npix = map_bin.npix |
|
101
|
|
|
see_pixels = map_bin.observation_map._pixels_ids |
|
102
|
|
|
transits = map_bin.n_transits |
|
103
|
|
|
scheme = map_bin.scheme |
|
104
|
|
|
|
|
105
|
|
|
nest_scheme=False |
|
106
|
|
|
if scheme.lower()=='nested': |
|
107
|
|
|
nest_scheme=True |
|
108
|
|
|
|
|
109
|
|
|
#what we want |
|
110
|
|
|
data = map_bin.observation_map.as_dense() |
|
111
|
|
|
bkg = map_bin.background_map.as_dense() |
|
112
|
|
|
|
|
113
|
|
|
zeros = np.empty(npix) |
|
114
|
|
|
zeros.fill(9e9) |
|
115
|
|
|
|
|
116
|
|
|
outFileName="{0}_bin{1}.fits.gz".format(outfile,analysis_bin) |
|
117
|
|
|
#I probably need to add some header info, but yeah |
|
118
|
|
|
hp.fitsfunc.write_map(outFileName, (data, bkg, zeros), |
|
119
|
|
|
column_names=labels, column_units=label_units, dtype=label_format, |
|
120
|
|
|
partial=False, fits_IDL=True, overwrite=clobber, nest=nest_scheme) |
|
121
|
|
|
|
|
122
|
|
|
|
|
123
|
|
|
#add the cards to the header |
|
124
|
|
|
with fits.open(outFileName,'update') as hdu1: |
|
125
|
|
|
hdr=hdu1[0].header |
|
126
|
|
|
|
|
127
|
|
|
for i,key in enumerate(primary_keys): |
|
128
|
|
|
val=primary_values[i] |
|
129
|
|
|
comment=primary_comments[i] |
|
130
|
|
|
|
|
131
|
|
|
if key=='TOTDUR': |
|
132
|
|
|
val=24.0*transits |
|
133
|
|
|
|
|
134
|
|
|
if key=='STOPMJD': |
|
135
|
|
|
val=startMJD+transits |
|
136
|
|
|
|
|
137
|
|
|
entry=(val,comment) |
|
138
|
|
|
|
|
139
|
|
|
hdr[key]=entry |
|
140
|
|
|
|
|
141
|
|
|
#File is now closed |
|
142
|
|
|
print("File Written: {0}".format(outFileName)) |
|
143
|
|
|
|
giacomov /
hawc_hal
