-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathplot_peaks.py
executable file
·107 lines (83 loc) · 3.33 KB
/
plot_peaks.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
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
#!/usr/bin/env python
"""
# 02_plot_spectra.py
Plots the time-averaged spectra.
"""
import matplotlib as mpl
import os
import seaborn as sns
import tables as tb
from leda_cal.leda_cal import *
from leda_cal.skymodel import *
from leda_cal.dpflgr import *
sns.set_style('ticks')
sns.set_context("paper",font_scale=1.5)
def quicklook(filename, save, flag, noshow):
h5 = tb.open_file(filename)
T_ant = apply_calibration(h5)
f_leda = T_ant['f']
ant_ids = ['252A', '254A', '255A', '252B', '254B', '255B'] # Added 252B, why was it missing? HG
mmax = 0
for ant in ant_ids:
if flag:
T_flagged = rfi_flag(T_ant[ant], thr_f=0.2, thr_t=0.2, rho=1.5,
bp_window_f=16, bp_window_t=16,
max_frac_f=0.5, max_frac_t=0.5)
T_ant[ant] = np.ma.filled(T_flagged, 0)
if flag: tmax = np.max(T_ant[ant])
else: tmax = np.percentile(T_ant[ant], 99.95)
if tmax > mmax: mmax = tmax
print("Plotting...")
fig, ax = plt.subplots(figsize=(8, 6))
mid = T_ant["252A"].shape[0]/2 # Currently plots the middle bit of the observation
print T_ant['lst'][mid] # Print the LST about which we are averaging
sl = 50 # The number of integrations to average either side
print ant_ids[0], mid, sl
print ant_ids[1], mid, sl
print ant_ids[2], mid, sl
plt.subplot(2,1,1)
plt.plot(f_leda, np.max(T_ant[ant_ids[0]], axis=0), label=ant_ids[0])
plt.plot(f_leda, np.max(T_ant[ant_ids[1]], axis=0), label=ant_ids[1])
plt.plot(f_leda, np.max(T_ant[ant_ids[2]], axis=0), label=ant_ids[2])
plt.legend(frameon=False)
plt.ylabel("Temperature [K]")
plt.xlim(30, 85)
plt.minorticks_on()
plt.ylim(500, mmax)
plt.subplot(2,1,2)
plt.plot(0, 0)
plt.plot(f_leda, np.max(T_ant[ant_ids[3]], axis=0), label=ant_ids[3])
plt.plot(f_leda, np.max(T_ant[ant_ids[4]], axis=0), label=ant_ids[4])
plt.plot(f_leda, np.max(T_ant[ant_ids[4]], axis=0), label=ant_ids[5])
plt.legend(frameon=False)
plt.xlim(30, 85)
plt.minorticks_on()
plt.ylim(500, mmax)
plt.xlabel("Frequency [MHz]")
plt.ylabel("Temperature [K]")
plt.legend(frameon=False)
plt.tight_layout()
plt.savefig("figures/compare-spectra.pdf")
if save:
if flag: plt.savefig("peaks_"+os.path.basename(filename)[:-3]+"_"+"flagged.png")
else: plt.savefig("peaks_"+os.path.basename(filename)[:-3]+".png")
if not noshow:
plt.show()
if __name__ == "__main__":
import optparse, sys
usage = '%prog [opts] filename_of_hdf5_observation'
o = optparse.OptionParser()
o.set_usage(usage)
o.set_description(__doc__)
o.add_option('--flag', dest='flag', action='store_true', default=False,
help='Apply flagging. Default: False')
o.add_option('--save', dest='save', action='store_true', default=False,
help="Save the plot to an image file, with filename the same as the h5 but png extension. Default: False.")
o.add_option('--noshow', dest='noshow', action='store_true', default=False,
help="Don't display the plot on screen. Useful for batch runs. Default: False.")
opts, args = o.parse_args(sys.argv[1:])
if len(args) != 1:
o.print_help()
exit(1)
else: filename = args[0]
quicklook(filename, opts.save, opts.flag, opts.noshow)