flatframes.py 12.7 KB
Newer Older
Monica Rainer's avatar
Monica Rainer committed
1
2
3
4
"""
Last modified: 2017-03-07

Reduction of the flat-fields:
Monica Rainer's avatar
Monica Rainer committed
5
- NOT ANYMORE: check the exposure time
Monica Rainer's avatar
Monica Rainer committed
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
- check the flux in a chosen area
- subtract the masterdark
- remove bad pixel using the bad pixel mask
- if more than one flats, combine them and adjust ron and gain
- save the masterflat in the calibration database
- straigthen the masterflat
- normalize the masterflat
- use the masterflat to create the masked array for extraction
"""

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

#import matplotlib.pyplot as plt


class GBFlats():
    def __init__(self, flats, dbconn):
        self.flats = flats
        self.dbconn = dbconn
        self.quality = []
        self.messages = []
        self.flatlist = []
        self.flatcorr = []


    def qualitycheck(self):
        """
        Check flat's quality: exposure time and signal in a well-defined region
        """

        for frame in self.flats:
            #print frame

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

# no exposure time check at the moment

#            if flat.header[CONFIG['KEYS']['EXPTIME']] != CONFIG['FLAT_EXPT']:
#                self.messages.append('Flat frame %s has unexpected exposure time of %s sec: skipped.' % (str(os.path.basename(frame)),str(flat.header[CONFIG['KEYS']['EXPTIME']]),))
#                continue


# check the signal in a well defined zone

            #zone = flat.data[CONFIG['FLATCHECK'][0]:CONFIG['FLATCHECK'][1],CONFIG['FLATCHECK'][2]:CONFIG['FLATCHECK'][3]]
            #mean = np.mean(zone)
            mean = np.mean(flat.data)


            self.messages.append('Mean value of flat-field %s: %s ' % (os.path.basename(frame),str(mean),))

            if mean < CONFIG['FLATSIGNAL']:
                self.messages.append('Flat frame %s failed quality check: signal too low (%s)' % (str(os.path.basename(frame)),str(mean)),)
                flat = None
69
                self.quality.append('FAILED')
Monica Rainer's avatar
Monica Rainer committed
70
71
            else:
                self.flatlist.append(flat)
72
                self.quality.append('OK')
Monica Rainer's avatar
Monica Rainer committed
73
74
75
76
77
78
79
80
        return

    def masterflat(self):
        """
        Create masterflat: remove dark, remove bad pixels, combine the images
        using ron and gain, correct ron and gain --> ron_eff and gain_eff
        Write the result in the database.
        """
81

Monica Rainer's avatar
Monica Rainer committed
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
        if self.flatlist:
            exptime = self.flatlist[0].header[CONFIG['KEYS']['EXPTIME']]
            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 flat-field. 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 flat in self.flatlist:

121
                # subdtract the masterdark
Monica Rainer's avatar
Monica Rainer committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
                if use_dark:
                    subflat = ccdproc.subtract_dark(flat,mdark,exposure_time=CONFIG['KEYS']['EXPTIME'],exposure_unit=u.second)
                    self.messages.append('Masterdark subtracted.')
                else:
                    subflat = flat

                # mask the flat using the badpix mask
                flat_corrected = varie.badpix(subflat.data,bad_mask,inverse_mask)
                self.messages.append('Bad pixel correction done.')

                flat_corr = ccdproc.CCDData(flat_corrected,unit=u.adu)

                self.flatcorr.append(flat_corr)

Monica Rainer's avatar
Monica Rainer committed
136
            if len(self.flatcorr) > 1:
Monica Rainer's avatar
Monica Rainer committed
137

Monica Rainer's avatar
Monica Rainer committed
138
139
140
141
142
143
144
145
                combineflat = ccdproc.Combiner(self.flatcorr)

                # mask the pixels using ron/gain
                #combineflat.sigma_clipping(func=np.ma.mean, dev_func=varie.stdcombine)
                combineflat.sigma_clipping(func=np.ma.median, dev_func=varie.stdcombine)
                mflat = combineflat.average_combine()
            else:
                mflat = self.flatcorr[0]
Monica Rainer's avatar
Monica Rainer committed
146

Monica Rainer's avatar
Monica Rainer committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
            mflat.data = np.asarray(mflat.data, dtype='float32')

            hea = self.flatlist[0].header
            flatname = hea[CONFIG['KEYS']['IMANAME']].replace('.fts','_FLAT.fits')

            nome = os.path.join(CONFIG['CALIB_DIR'],flatname)

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

            hdu = fits.PrimaryHDU(data=mflat.data,header=mflat.header)
            masterflat = fits.HDUList([hdu])
            masterflat.writeto(nome,clobber=True)

            db.insert_dbfile(self.dbconn,'flat',nome)
            self.messages.append('Masterflat %s was created and inserted in the calibration database' % (str(os.path.basename(nome))))

167
168
169
170
171
172
            shiftY = varie.shiftY(mflat.data)

            db.insert_dbfile(self.dbconn,'shiftY',shiftY)

            mflat = None

Monica Rainer's avatar
Monica Rainer committed
173
174
175
176
177
        else:
            self.messages.append('There are no available flats. The masterflat will be taken from the calibration database.')
            if db.check_dbfile(self.dbconn,'flat'):
                try:
                    db.copy_dbfile(self.dbconn,'flat')
178
179
180
181
182
183
184
185
186
187
188
189
190
191
                    try:
                        db.copy_dbfile(self.dbconn,'shiftY')
                        shift = db.extract_dbfile(self.dbconn,'shiftY')
                        if not shift:
                            cal_flat = db.extract_dbfile(self.dbconn,'flat')
                            mflat = ccdproc.CCDData.read(cal_flat, unit=u.adu)
                            shiftY = varie.shiftY(mflat.data)
                            db.insert_dbfile(self.dbconn,'shiftY',shiftY)

                    except:
                        cal_flat = db.extract_dbfile(self.dbconn,'flat')
                        mflat = ccdproc.CCDData.read(cal_flat, unit=u.adu)
                        shiftY = varie.shiftY(mflat.data)
                        db.insert_dbfile(self.dbconn,'shiftY',shiftY)
Monica Rainer's avatar
Monica Rainer committed
192
193
194
                except:
                    self.messages.append('There are no flats in the calibration database.')
                    return False
195
196
197
198
199
200
201

            if db.check_dbfile(self.dbconn,'shiftY'):
                cal_flat = db.extract_dbfile(self.dbconn,'flat')
                mflat = ccdproc.CCDData.read(cal_flat, unit=u.adu)
                shiftY = varie.shiftY(mflat.data)
                db.insert_dbfile(self.dbconn,'shiftY',shiftY)

Monica Rainer's avatar
Monica Rainer committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
        return True


    def straighten(self):
        """
        Straighten the masterflat, normalize it and create the masked array
        for A, B and C extraction
        """
        # straighten the flat

        flat = db.extract_dbfile(self.dbconn,'flat')

        straight = flat.replace('.fits','_str.fits')
        nor = flat.replace('.fits','_nor.fits')
Monica Rainer's avatar
Monica Rainer committed
216
        #extracted = flat.replace('.fits','_ext.fits')
Monica Rainer's avatar
Monica Rainer committed
217

218
219
220
221
222
223
224
225
226
227
228
229
        args = [CONFIG['STRAIGHT'],flat,straight]
        args.extend(CONFIG['STRAIGHT_OPT'])
        # search for shift defined in the straighten options in config.py
        dy = True
        for opt in CONFIG['STRAIGHT_OPT']:
            try:
                dy = opt.rindex('DY=')
                ypos = int(opt[dy-2:])
                shift = CONFIG['SHIFT_Y'] + ypos
                dy = False
            except:
                pass
230

231
232
233
234
235
236
237
        if dy:
            shift = db.extract_dbfile(self.dbconn,'shiftY')
            if not shift:
                cal_flat = db.extract_dbfile(self.dbconn,'flat')
                mflat = ccdproc.CCDData.read(cal_flat, unit=u.adu)
                shift = varie.shiftY(mflat.data)
                db.insert_dbfile(self.dbconn,'shiftY',shift)
238

239
240
241
            shiftY = [''.join(('DY=',str(shift - CONFIG['SHIFT_Y'])))]
            #print shiftY
            args.extend(shiftY)
242

243
        #print args
Monica Rainer's avatar
Monica Rainer committed
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
        subprocess.call(args)

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

        self.messages.append('Orders straightened.')
        db.insert_dbfile(self.dbconn,'flatstr',straight)

        mflat = ccdproc.CCDData.read(straight, unit=u.adu)
        fdata = mflat.data

        norflat = np.ones((CONFIG['YCCD'],CONFIG['XCCD']), dtype='float32')
        #maskA = np.ones((CONFIG['YCCD'],CONFIG['XCCD']), dtype=np.int)
        #maskB = np.ones((CONFIG['YCCD'],CONFIG['XCCD']), dtype=np.int)
        maskC = np.ones((CONFIG['YCCD'],CONFIG['XCCD']), dtype='int')

        xvalues = np.arange(CONFIG['XCCD'])
        np.seterr(divide='ignore')

Monica Rainer's avatar
Monica Rainer committed
263
264
        #ordini = np.zeros((CONFIG['N_ORD'],CONFIG['XCCD']))

Monica Rainer's avatar
Monica Rainer committed
265
        for x in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
266
267
268
269

            #frow = x*CONFIG['W_ORD'] + 19
            #ordini[x] = np.sum(fdata[frow-5:frow+5],axis=0)

Monica Rainer's avatar
Monica Rainer committed
270
            # order limits
Monica Rainer's avatar
Monica Rainer committed
271
272
            start = x*CONFIG['W_ORD']
            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
Monica Rainer's avatar
Monica Rainer committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298

            # evaluate average background value for each order
            #background = []
            #background.append(np.median(fdata[start+2]))
            #background.append(np.median(fdata[start+3]))
            #background.append(np.median(fdata[end-2]))
            background = np.sort(fdata,axis=None)[0:CONFIG['YCCD']*3]
            back = np.median(background)

            # create mask for extraction

            #i = 0
            #width = 0
            for row in xrange(3,38,1):
                rrow = start + row
                if np.median(fdata[rrow]) > back*2:
                #if np.median(fdata[rrow]) > back*1.7:
                    fit = np.polyfit(xvalues,fdata[rrow],5)
                    pfit = np.polyval(fit,xvalues)
                    norflat[rrow] = fdata[rrow]/pfit
                    maskC[rrow] = 0
                    #if not i:
                    #    i = rrow
                    #width = width + 1


Monica Rainer's avatar
Monica Rainer committed
299
300
301
302
303
304
305
306
307
308
        # save the extracted flat as FITS file
        #flextract = ccdproc.CCDData(ordini, unit=u.adu)
        #flextract.data = np.asarray(flextract.data, dtype='float32')
        #flextract.header = mflat.header
        #fhdu = fits.PrimaryHDU(data=flextract.data,header=flextract.header)
        #flextr = fits.HDUList([fhdu])
        #flextr.writeto(extracted,clobber=True)
        #red_ext = os.path.join(CONFIG['RED_CALIB'],os.path.basename(extracted))
        #shutil.copyfile(extracted,red_ext)

Monica Rainer's avatar
Monica Rainer committed
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
        fnor = ccdproc.CCDData(norflat, unit=u.adu)
        fnor.header = mflat.header
        #fn = fnor.to_hdu()
        hdu = fits.PrimaryHDU(data=fnor.data,header=fnor.header)
        fn = fits.HDUList([hdu])
        fn.writeto(nor,clobber=True)

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

        db.insert_dbfile(self.dbconn,'flatnor',nor)
        self.messages.append('Flat-field normalized.')
        self.messages.append('Normalized masterflat %s was created and inserted in the calibration database' % (str(os.path.basename(nor))))

        #cmask = ccdproc.CCDData(maskC, unit=u.adu)
        cmask = np.asarray(maskC, dtype='int')
        #cm = cmask.to_hdu()
        hdu = fits.PrimaryHDU(data=cmask)
        cm = fits.HDUList([hdu])
        cm.writeto(CONFIG['MASK_C'],clobber=True)

        self.messages.append('The extraction mask was created.')
        return


    def process(self):
        warnings.simplefilter('ignore', category=AstropyWarning)
        self.qualitycheck()
        if self.masterflat():
            self.straighten()
        #self.flats[:] = []
        return