wlframes.py 10.7 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
"""
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)
49
            self.quality.append('OK')
Monica Rainer's avatar
Monica Rainer committed
50
51
52
53
54
        return

    def masterlamp(self,unefp):

        if self.wllist:
55
            exptime = float(self.wllist[0].header[CONFIG['KEYS']['EXPTIME']])
Monica Rainer's avatar
Monica Rainer committed
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
            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
109
110
            #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
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

            #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')

160
161
162
        #args = [CONFIG['STRAIGHT'],lamp,straight,CONFIG['STRAIGHT_OPT']]
        args = [CONFIG['STRAIGHT'],lamp,straight]
        args.extend(CONFIG['STRAIGHT_OPT'])
Monica Rainer's avatar
Monica Rainer committed
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
        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
179
        ordini = np.zeros((CONFIG['N_ORD'],CONFIG['XCCD']))
Monica Rainer's avatar
Monica Rainer committed
180
181
182
183
184
185
186

        # 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
187
188
        for x in xrange(CONFIG['N_ORD']):
            row = x*CONFIG['W_ORD'] + 19
Monica Rainer's avatar
Monica Rainer committed
189
 
Monica Rainer's avatar
Monica Rainer committed
190
            ordini[x] = np.sum(wldata[row-CONFIG['WEXT']:row+CONFIG['WEXT']],axis=0)
Monica Rainer's avatar
Monica Rainer committed
191
            if unefp == 'une':
192
                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
193
194
195
196
197
198
199
                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:
200
                        wcalib = db.extract_dbfile(self.dbconn,'une_calib')
Monica Rainer's avatar
Monica Rainer committed
201
202
203
204
205
206
                    except:
                        wcalib = False

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

207
208
                        db.copy_dbfile(self.dbconn,'une_calib')
                        wcalib = db.extract_dbfile(self.dbconn,'une_calib')
Monica Rainer's avatar
Monica Rainer committed
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224

                    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:
Monica Rainer's avatar
Monica Rainer committed
225
                    wcalib = db.extract_dbfile(self.dbconn,'une_calib')
Monica Rainer's avatar
Monica Rainer committed
226
227
228
229
230
231
                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:
Monica Rainer's avatar
Monica Rainer committed
232
233
                        db.copy_dbfile(self.dbconn,'une_calib')
                        wcalib = db.extract_dbfile(self.dbconn,'une_calib')
Monica Rainer's avatar
Monica Rainer committed
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
278
279
280
                    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