-
Notifications
You must be signed in to change notification settings - Fork 3
/
inspectSED.py
executable file
·213 lines (188 loc) · 7.58 KB
/
inspectSED.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
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
209
210
211
212
#!/usr/bin/env python3
import argparse
from sed_input import read_ascii, read_cleaned, read_spectrum, magToJy, JyToLamFlam
import sys, os
import matplotlib.pyplot as plt
from matplotlib.pyplot import errorbar, loglog
import numpy as np
from plot import pltSED
from pathlib import Path
description = \
"""
description:
Read in photometric and spectroscopic data and plot it
so that spurious data points may be rejected
from the final data file.
"""
epilog = \
"""
examples:
inspectPhot.py --phot=AKSco/AKSco_phot.dat
--spec=AKSco/28902101_sws.fit,AKSco/cassis_yaaar_spcfw_12700160t.fits
--pltR=0.1,1000
"""
parser = argparse.ArgumentParser(description=description,epilog=epilog,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("--phot",dest="phot",default='',type=str,
help='Full path to photometry data file.')
parser.add_argument("--spec",dest='spec',default='',type=str,
help='Full path to the Spitzer spectrum file.')
parser.add_argument("--scale",dest='specScale',default='',type=str,
help='Scale factor to be applied to spectral flux')
parser.add_argument("--pltR",dest='plt_range',default='[]',type=str,
help='X-range (in microns) for plot window (free by default)')
parser.add_argument("--savePlt",dest='saveplt',default=False,type=bool,
help='Save a .pdf copy of the full and cleaned SEDs (default False)')
argopt = parser.parse_args()
if argopt.specScale != '':
if len(argopt.specScale.split(',')) != len(argopt.spec.split(',')):
print('Error: options parsed to --spec and to --scale ')
print('must have the same length!')
sys.exit()
############
# 1. Read in the photometric data:
############
infile = Path(argopt.phot)
if infile.suffix == '.dat':
if 'cleaned' not in infile.name:
wvlen,wband,jy,ejy,flag,unit,beam,odate,ref = read_ascii(infile)
else:
wvlen,wband,f,ef,flag,beam,odate,ref = read_cleaned(infile)
jy = None
else:
print('')
print('File name error: function limited to plotting ascii files output by queryDB.py.')
print('')
sys.exit()
# If provided, read in the spectroscopic data:
if argopt.spec != '':
specFiles = argopt.spec.split(',')
if argopt.specScale != '':
specS = [float(s) for s in argopt.specScale.split(',')]
else:
specS = [1]*len(specFiles)
else:
specFiles = None
specS = None
############
# 2. Convert photometry data to W/m^2 (lamFlam):
############
if jy:
for i in range(0, len(jy)):
if ejy[i] == '--':
ejy[i] = np.nan
else:
ejy[i] = float(ejy[i])
if unit[i] == 'mag':
jy[i], ejy[i] = magToJy(jy[i],ejy[i], wband[i])
unit[i] = 'Jy'
elif unit[i] == 'mJy':
jy[i] = jy[i]*1e-3
ejy[i] = ejy[i]*1e-3
unit[i] = 'Jy'
f, ef = [], []
for i in range(0, len(wvlen)):
f.append(JyToLamFlam(jy[i],ejy[i],wvlen[i])[0])
ef.append(JyToLamFlam(jy[i],ejy[i],wvlen[i])[1])
############
# 5. Plot SED:
############
if argopt.plt_range != '[]':
try:
x_range = argopt.plt_range.split(',')
except:
print('Error: format of plot limits not recognised')
print('Defaults will be used instead')
x_range = 'default'
else:
x_range = 'default'
if argopt.saveplt == True:
pltSED(infile, x_range, f, ef, wvlen, specFiles, specS, interactive=False)
sedOutF = infile.parent / Path(infile.name.split('_')[0]+'_sed.pdf')
k = 0
while (sedOutF.parent / Path(sedOutF.name.replace('sed.pdf', 'sed_'+str(k)+'.pdf'))).exists():
k += 1 # avoids over-writing existing files
sedOutF = sedOutF.parent / Path(sedOutF.name.replace('sed.pdf', 'sed_'+str(k)+'.pdf'))
plt.savefig(sedOutF)
print('')
print('-----------------------------------------------')
print('| The plot displays your input SED data. |')
print('| |')
if argopt.saveplt == True:
print('| A copy of this plot has been saved to |')
print('| '+str(sedOutF))
print('| |')
print('| Please click on any photometry points you |')
print('| wish to remove from the final data set. |')
print('| |')
print('| The x and y positions, together with the |')
print('| waveband of each data point you click will |')
print('| be printed to the terminal screen. |')
print('| |')
print('| When you are finished, please close the |')
print('| plot window. |')
print('-----------------------------------------------')
fig = pltSED(infile, x_range, f, ef, wvlen, specFiles, specS, interactive=True)
indices = []
def on_pick(event):
global indices
artist = event.artist
xmouse, ymouse = event.mouseevent.xdata, event.mouseevent.ydata
x, y = artist.get_xdata(), artist.get_ydata()
ind = event.ind
indices.append(ind[0])
print('')
print('Detected mouse click at:')
print('x=', x[ind[0]], 'm; y=', y[ind[0]], 'W/m^2')
print('Corresponding to waveband:',wband[ind[0]])
print('')
plt.connect('pick_event', on_pick)
plt.show()
############
# 5. Write retained photometric data to file:
############
if 'cleaned' not in infile.name or indices != []:
print('')
print('Writing cleaned data to new file:')
if 'cleaned' not in infile.name:
outfile = infile.parent / Path(infile.stem+'_cleaned_')
else:
outfile = infile.parent / Path(infile.stem)
j = 0
while Path(str(outfile)+str(j)+'.dat').exists():
j += 1 # avoids over-writing other attempts to clean data file
outfile = Path(str(outfile)+str(j)+'.dat')
print(outfile) # name of file to be written
print('')
with open(outfile,'w') as f_out:
with open(infile, 'r') as f_in:
for line in f_in:
try:
a = float(line[0])
except ValueError:
f_out.write(line.replace('mag','lamFlam').replace('m -- -- --','m -- W/m^2 W/m^2'))
for i in range(0, len(wvlen)):
if i not in indices:
f_out.write(' '.join([str(x) for x in [wvlen[i], # wavelength in m
wband[i], # waveband
f[i], # flux in W/m^2
ef[i], # flux error in W/m^2
flag[i], # flag on flux
beam[i], # beam size (may be dummy value)
odate[i], # obs date (may be unknown or average)
ref[i]]])+'\n') # reference for original data
############
# 6. Save cleaned version of SED plot to file
############
if argopt.saveplt == True and indices != []:
# read in cleaned data:
wvlen,wband,f,ef,flag,beam,odate,ref = read_cleaned(outfile)
pltSED(infile, x_range, f, ef, wvlen, specFiles, specS, interactive=False)
h = 0
while (sedOutF.parent / Path(sedOutF.name.replace('sed_'+str(k)+'.pdf', 'sed_cleaned_'+str(h)+'.pdf'))).exists():
h += 1 # avoids over-writing existing files
sedOutF = sedOutF.parent / Path(sedOutF.name.replace('sed_'+str(k)+'.pdf', 'sed_cleaned_'+str(h)+'.pdf'))
plt.savefig(sedOutF)
print('')
print('Info: Plot of cleaned SED saved to '+str(sedOutF))
print('')