wlframes.py 10.6 KB
Newer Older
Monica Rainer's avatar
Monica Rainer committed
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
"""
Last modified: 2017-03-07

Reduction of the UNe lamps:
- check the exposure time???
- check the flux ratio in chosen areas???
- subtract the masterdark
- remove bad pixel using the bad pixel mask
- if more than one lamps, combine them and adjust ron and gain
- save the masterlamp in the calibration database
- straigthen the masterlamp
- extract the central 4 pixels of each order for quality check
- UNe wavelength calibration (TO DO: FP)
- TO DO: UNe line diagnostics (checking the ratio between specific emission lines)
"""

from astropy import units as u
from astropy.io import fits
import warnings
from astropy.utils.exceptions import AstropyWarning


import ccdproc
from drslib.config import CONFIG
from drslib import db, varie
import numpy as np
import math, os, subprocess, shutil


class GBWls():
    def __init__(self, wls, dbconn):
        self.wls = wls
        self.dbconn = dbconn
        self.quality = []
        self.messages = []
        self.wllist = []
        self.wlcorr = []


    def qualitycheck(self,unefp):
        """
        Quality check for UNe and FP lamps still to be defined.
        """
        for frame in self.wls:

            wl = ccdproc.CCDData.read(frame, unit=u.adu)

            self.wllist.append(wl)
        return

    def masterlamp(self,unefp):

        if self.wllist:
54
            exptime = float(self.wllist[0].header[CONFIG['KEYS']['EXPTIME']])
Monica Rainer's avatar
Monica Rainer committed
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
            darkname = 'dark' + str(int(exptime))
            use_dark = True

            try:
                masterdark = db.extract_dbfile(self.dbconn,darkname)
            except:
                masterdark = False

            if not masterdark:
                self.messages.append('No masterdark found for this night, it will be taken from the calibration database.')

                try:
                    db.copy_dbfile(self.dbconn,darkname)
                    masterdark = db.extract_dbfile(self.dbconn,darkname)
                except:
                    for key in CONFIG['DARKLIST']:
                        try:
                            darkname = 'dark' + str(int(key))
                            db.copy_dbfile(self.dbconn,darkname)
                            masterdark = db.extract_dbfile(self.dbconn,darkname)
                            self.messages.append('No masterdark in the calibration database with the same exposure time as the lamp. The %s sec masterdark will be used instead' % (str(int(key))))
                            break
                        except:
                            self.messages.append('There are no masterdark in the calibration database. The masterdark will not be used.')
                            use_dark = False

            if use_dark:
                mdark = ccdproc.CCDData.read(masterdark, unit=u.adu)
            badpix = ccdproc.CCDData.read(CONFIG['BADPIX_MASK'], unit=u.adu)
            bad_mask=badpix.data
            inverse_mask=np.logical_not(bad_mask)


            for wl in self.wllist:

                # subdtract the masterdak
                if use_dark:
                    subwl = ccdproc.subtract_dark(wl,mdark,exposure_time=CONFIG['KEYS']['EXPTIME'],exposure_unit=u.second)
                    self.messages.append('Masterdark subtracted.')
                else:
                    subwl = wl

                # mask the lamp using the badpix mask
                wl_corrected = varie.badpix(subwl.data,bad_mask,inverse_mask)
                self.messages.append('Bad pixel correction done.')

                wl_corr = ccdproc.CCDData(wl_corrected,unit=u.adu)

                self.wlcorr.append(wl_corr)

            combinewl = ccdproc.Combiner(self.wlcorr)

            # mask the pixels using ron/gain
108
109
            #combinewl.sigma_clipping(func=np.ma.mean, dev_func=varie.stdcombine)
            combinewl.sigma_clipping(func=np.ma.median, dev_func=varie.stdcombine)
Monica Rainer's avatar
Monica Rainer committed
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

            #combineflat.sigma_clipping(func=np.ma.mean)
            mwl = combinewl.average_combine()
            mwl.data = np.asarray(mwl.data, dtype='float32')

            hea = self.wllist[0].header
            if unefp == 'une':
                wlname = hea[CONFIG['KEYS']['IMANAME']].replace('.fts','_UNE.fits')

            else:
                wlname = hea[CONFIG['KEYS']['IMANAME']].replace('.fts','_FP.fits')

            wlname = wlname.replace(' ', '').replace('-', '_')
            nome = os.path.join(CONFIG['CALIB_DIR'],wlname)

            mwl.header = hea
            mwl.header[CONFIG['KEYS']['FILENAME']] = wlname
            mwl.header[CONFIG['KEYS']['NCOMBINE']] = len(self.wllist)
            mwl.header[CONFIG['RON_EFF'][0]] = (math.sqrt(len(self.wllist))*CONFIG['RON'],CONFIG['RON_EFF'][1])
            mwl.header[CONFIG['GAIN_EFF'][0]] = (len(self.wllist)*CONFIG['GAIN'],CONFIG['GAIN_EFF'][1])

            hdu = fits.PrimaryHDU(data=mwl.data,header=mwl.header)
            masterwl = fits.HDUList([hdu])
            masterwl.writeto(nome,clobber=True)
            mwl = None
            db.insert_dbfile(self.dbconn,unefp,nome)
            self.messages.append('%s masterlamp %s was created and inserted in the calibration database.' % (unefp,str(os.path.basename(nome))))


        else:
            self.messages.append('There are no available %s lamps. The reference lamp will be taken from the calibration database.' % (unefp))

            if db.check_dbfile(self.dbconn,unefp):
                try:
                    db.copy_dbfile(self.dbconn,unefp)
                except:
                    self.messages.append('There are no %s lamps in the calibration database.' %(unefp))
                    return False
        return True


    def extract(self,unefp):
        # straighten the lamp's orders
        lamp = db.extract_dbfile(self.dbconn,unefp)

        #print lamp
        straight = lamp.replace('.fits','_str.fits')
        extracted = lamp.replace('.fits','_ext.fits')

        args = [CONFIG['STRAIGHT'],lamp,straight,CONFIG['STRAIGHT_OPT']]
        subprocess.call(args)

        red_straight = os.path.join(CONFIG['RED_CALIB'],os.path.basename(straight))
        shutil.copyfile(straight,red_straight)

        #self.messages.append('%s: orders straightened.' % (str(os.path.basename(lamp))))

        strname = '_'.join((unefp,'str'))
        db.insert_dbfile(self.dbconn,strname,straight)
        self.messages.append('Straightened %s masterlamp %s was created and inserted in the calibration database.' % (unefp,str(os.path.basename(straight))))

        wl = ccdproc.CCDData.read(straight, unit=u.adu)

        wldata = wl.data
        wl.header[CONFIG['WLFIT'][0]] = (CONFIG['WLFIT_FUNC'],CONFIG['WLFIT'][1])

Monica Rainer's avatar
Monica Rainer committed
176
        ordini = np.zeros((CONFIG['N_ORD'],CONFIG['XCCD']))
Monica Rainer's avatar
Monica Rainer committed
177
178
179
180
181
182
183

        # read the lines to use in the wavelength calibration
        select_lines, all_lines = varie.UNe_linelist() 

        # extract the central 10 pixel of each order
        # wavelength calibration, to check the lamp and
        # to be used as default if later ad hoc calibration fails
Monica Rainer's avatar
Monica Rainer committed
184
185
        for x in xrange(CONFIG['N_ORD']):
            row = x*CONFIG['W_ORD'] + 19
Monica Rainer's avatar
Monica Rainer committed
186
187
188
 
            ordini[x] = np.sum(wldata[row-5:row+5],axis=0)
            if unefp == 'une':
189
                calib_failed, coeffs, comments = varie.UNe_calibrate(ordini[x],x+32,select_lines[x+32],all_lines[x+32],True,True)
Monica Rainer's avatar
Monica Rainer committed
190
191
192
193
194
195
196
                for comment in comments:
                    self.messages.append(comment)

                if calib_failed:
                    self.messages.append(' *** WARNING ***')
                    self.messages.append('The default wavelength calibration for the order %s will be taken from the database and as such it will not be optimal for the night.' % str(x+32))
                    try:
197
                        wcalib = db.extract_dbfile(self.dbconn,'une_calib')
Monica Rainer's avatar
Monica Rainer committed
198
199
200
201
202
203
                    except:
                        wcalib = False

                    if not wcalib:
                        self.messages.append('No UNe calibration found for this night, it will be taken from the calibration database.')

204
205
                        db.copy_dbfile(self.dbconn,'une_calib')
                        wcalib = db.extract_dbfile(self.dbconn,'une_calib')
Monica Rainer's avatar
Monica Rainer committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277

                    try:
                        wlc = ccdproc.CCDData.read(wcalib, unit=u.adu)
                        wl.header = wlc.header
                    except:
                        self.messages.append(' *** WARNING ***')
                        self.messages.append('There is no default wavelength calibration in the database, please observe a new lamp, otherwise the spectra will not be wavelength calibrated.')

                else:
                    for key in coeffs:
                        keyword = ''.join((CONFIG['WLCOEFFS'][key][0],str(x+32)))
                        wl.header[keyword] = (coeffs[key],CONFIG['WLCOEFFS'][key][1])

            else:
                # calibrate FP using UNe
                try:
                    wcalib = db.extract_dbfile(self.dbconn,'une_str')
                except:
                    wcalib = False

                if not wcalib:
                    self.messages.append('No UNe calibration found for this night, it will be taken from the calibration database.')
                    try:
                        db.copy_dbfile(self.dbconn,'une_str')
                        wcalib = db.extract_dbfile(self.dbconn,'une_str')
                    except:
                        return

                #print wcalib
                wlc = ccdproc.CCDData.read(wcalib, unit=u.adu)
                wl.header = wlc.header


        # save the extracted lamp as FITS file
        wlextract = ccdproc.CCDData(ordini, unit=u.adu)
        wlextract.data = np.asarray(wlextract.data, dtype='float32')
        wlextract.header = wl.header
        hdu = fits.PrimaryHDU(data=wlextract.data,header=wlextract.header)
        wlextr = fits.HDUList([hdu])
        wlextr.writeto(extracted,clobber=True)

        red_extr = os.path.join(CONFIG['RED_CALIB'],os.path.basename(extracted))
        shutil.copyfile(extracted,red_extr)

        calname = '_'.join((unefp,'calib'))
        db.insert_dbfile(self.dbconn,calname,extracted)
        self.messages.append('%s wavelength calibration %s was created and inserted in the calibration database.' % (unefp,str(os.path.basename(extracted))))

        wlextract = None
        ordini = None

        return


    def une_process(self):
        warnings.simplefilter('ignore', category=AstropyWarning)
        self.qualitycheck('une')
        if self.masterlamp('une'):
            self.extract('une')
        #self.wls[:] = []
        return

    def fp_process(self):
        warnings.simplefilter('ignore', category=AstropyWarning)
        self.qualitycheck('fp')
        if self.masterlamp('fp'):
            self.extract('fp')
        #self.wls[:] = []
        return