/
xsh_combine.py
119 lines (97 loc) · 3.87 KB
/
xsh_combine.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
108
109
110
111
112
113
114
115
116
117
118
119
# ------------------------------------------------------------------------------
# XSH_COMBINE
# Combine X-shooter spectra
# v2.0 - 2018-12-06
# Guido Cupani - INAF-OATs
# ------------------------------------------------------------------------------
# Sample run:
# > python xsh_combine.py -h // Help
# > python $APP_PATH/xsh_combine.py -l=xsh_combine_list.dat
# ------------------------------------------------------------------------------
import argparse
from astropy.io import ascii
from formats import XshMerge, save
from matplotlib import pyplot as plt
from matplotlib.figure import Figure
import numpy as np
from toolbox import air_to_vacuum, earth_to_bary
class ArmStack(object):
pass
def ave(data, wgh):
if len(np.shape(data)) == 1:
return data
if len(np.shape(data)) == 2:
return np.average(data, axis=0, weights=wgh)
def run(**kwargs):
# Load parameters
frames = np.array(ascii.read(kwargs['framelist'],
format='no_header')['col1'])
name = kwargs['name']
first = True
fig = Figure()
ax = plt.gca()
wave_stack = ArmStack()
flux_stack = ArmStack()
err_stack = ArmStack()
for f in frames:
print("Processing", f+"...")
dim = f[-11:-9]
arm = f[-8:-5]
# Load spectrum
spec = XshMerge(f)
# Correct wavelengths from air to vacuum
air_to_vacuum(spec)
spec.save(f[:-5]+'_VAC.fits', dim)
# Correct wavelengths from Earth to barycentric frame
earth_to_bary(spec)
spec.save(f[:-5]+'_VAC_BARY.fits', dim)
if dim == '1D':
#spec.err[spec.err == 0] = -9999.*1e-15
if hasattr(flux_stack, arm):
setattr(flux_stack, arm,
np.vstack((getattr(flux_stack, arm), spec.flux)))
setattr(err_stack, arm,
np.vstack((getattr(err_stack, arm), spec.err)))
else:
setattr(wave_stack, arm, spec.wave)
setattr(flux_stack, arm, spec.flux)
setattr(err_stack, arm, spec.err)
flux_ave = ArmStack()
err_ave = ArmStack()
spec_app = None
for arm in ['UVB', 'VIS', 'NIR']:
if hasattr(flux_stack, arm):
wave = getattr(wave_stack, arm)
flux = getattr(flux_stack, arm)
err = getattr(err_stack, arm)
setattr(flux_ave, arm, ave(flux, wgh=1/err**2))
setattr(err_ave, arm,
np.sqrt(ave(err**2, wgh=1/err**2)/np.shape(err)[0]))
spec_arr = np.array([wave, getattr(flux_ave, arm),
getattr(err_ave, arm)])
save(name+'_%s_VAC_BARY_COMB.fits' % arm, spec.hdr, wave,
getattr(flux_ave, arm), getattr(err_ave, arm))
if spec_app is not None:
spec_app = np.append(spec_app, spec_arr, axis=1)
else:
spec_app = spec_arr
ax.plot(spec_app[0], spec_app[1], c='black', label='Weighted average')
ax.plot(wave_stack.UVB, flux_ave.UVB+err_ave.UVB, c='black', linewidth=0.5)
ax.plot(wave_stack.UVB, flux_ave.UVB-err_ave.UVB, c='black', linewidth=0.5)
ax.legend()
plt.show()
def main():
""" Read the CL parameters and run """
p = argparse.ArgumentParser()
p.add_argument('-l', '--framelist', type=str,
default='xsh_combine_list.dat',
help="List of frames; must be an ascii file with a column "
"of entries. It's mandatory that frames are grouped in "
"UVB/VIS/NIR sets (in this order) and that at least 2 "
"groups are present.")
p.add_argument('-n', '--name', type=str,
default='noname', help="Name prefix for the output frames")
args = vars(p.parse_args())
run(**args)
if __name__ == '__main__':
main()