show_viscube.py
4.92 KB
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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
#! /usr/bin/env python
# ==========================================================================
# Show visibility cube
#
# Copyright (C) 2019 Juergen Knoedlseder
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# ==========================================================================
import sys
import gammalib
import cscripts
try:
import matplotlib.pyplot as plt
plt.figure()
plt.close()
except (ImportError, RuntimeError):
print('This script needs the "matplotlib" module')
sys.exit()
# =========================== #
# Convert map into Aitoff map #
# =========================== #
def get_aitoff_map(filename):
"""
Convert sky map file into an Aitoff map
Parameters
----------
filename : str
File name
Returns
-------
aitoff = `~gammalib.GSkyMap`
Aitoff sky map
"""
# Read sky map
map = gammalib.GSkyMap(filename)
# Integrate map over all zenith angles
map.stack_maps()
# Allocate Aitoff map
aitoff = gammalib.GSkyMap('AIT', 'GAL', 0.0, 0.0, -1.0, 1.0, 324, 162)
# Loop over all pixels in Aitoff map
for i in range(aitoff.npix()):
# Get sky direction
try:
# Get direction
dir = aitoff.inx2dir(i)
# Get value
value = map(dir)
# Set pixel value
aitoff[i] = value
except:
continue
# Return
return aitoff
# ==================== #
# Plot visibility cube #
# ==================== #
def plot_viscube(ax, filename):
"""
Plot visibility cube
Parameters
----------
ax : plt
Subplot
filename : str
File name
"""
# Get Aitoff map
map = get_aitoff_map(filename)
# Create array from skymap
array = []
for iy in range(map.ny()):
row = []
for ix in range(map.nx()):
index = ix+iy*map.nx()
value = map[index]
row.append(value)
array.append(row)
# Show Aitoff projection
c = ax.imshow(array,
extent=(gammalib.pi,-gammalib.pi,gammalib.pi/2,-gammalib.pi/2),
aspect=0.5)
cbar = plt.colorbar(c, ax=ax, orientation='vertical', shrink=0.7)
cbar.set_label('Visibility (hours)')
ax.grid(True)
ax.set_xlabel('Galactic longitude (deg)')
ax.set_ylabel('Galactic latitude (deg)')
ax.set_xticklabels([150, 120, 90, 60, 30, 0, 330, 300, 270, 240, 210])
# Return
return
# ======================= #
# Plot visibility results #
# ======================= #
def plot_visibility(ax, filename):
"""
Plot visibility results
Parameters
----------
ax : plt
Subplot
filename : str
File name
"""
# Load results
fits = gammalib.GFits(filename)
table = fits.table('VISIBILITY')
mjd = table['MJD']
dark_time = table['Darktime']
# Create arrays
nrows = table.nrows()
time = []
dark = []
for row in range(nrows):
time.append(mjd.real(row))
dark.append(dark_time.real(row))
# Plot arrays
ax.plot(time, dark, 'r-')
ax.set_xlabel('MJD (days)')
ax.set_ylabel('Dark time (hours)')
# Return
return
# ==================== #
# Show visibility cube #
# ==================== #
def show_viscube():
"""
Show visibility cube
"""
# Set usage string
usage = 'show_viscube.py [-t title] [-p plotfile] [filename]'
# Set default options
options = [{'option': '-p', 'value': ''},
{'option': '-t', 'value': 'Visibility map'}]
# Get arguments and options from command line arguments
args, options = cscripts.ioutils.get_args_options(options, usage)
# Extract script parameters from options
plotfile = options[0]['value']
title = options[1]['value']
# Create figure and subplots
fig = plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(211, projection='aitoff')
ax2 = fig.add_subplot(212)
fig.suptitle(title, fontsize=16)
# Plot visibility cube
plot_viscube(ax1, args[0])
# Plot visibility result
plot_visibility(ax2, args[0])
# Show figure
if len(plotfile) > 0:
plt.savefig(plotfile)
else:
plt.show()
# Return
return
# ======================== #
# Main routine entry point #
# ======================== #
if __name__ == '__main__':
# Show visibility cube
show_viscube()