flatframes.py 14.8 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
- 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
Andrea Bignamini's avatar
Andrea Bignamini committed
24
from drslib import metadata
Monica Rainer's avatar
Monica Rainer committed
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
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
70
                self.quality.append('FAILED')
Monica Rainer's avatar
Monica Rainer committed
71
72
            else:
                self.flatlist.append(flat)
73
                self.quality.append('OK')
Monica Rainer's avatar
Monica Rainer committed
74
75
76
77
78
79
80
81
        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.
        """
82

Monica Rainer's avatar
Monica Rainer committed
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
        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:

122
                # subdtract the masterdark
Monica Rainer's avatar
Monica Rainer committed
123
124
125
126
127
128
129
130
131
132
133
134
135
136
                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
137
            if len(self.flatcorr) > 1:
Monica Rainer's avatar
Monica Rainer committed
138

Monica Rainer's avatar
Monica Rainer committed
139
140
141
142
143
144
145
146
                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
147

Monica Rainer's avatar
Monica Rainer committed
148
149
150
151
152
153
154
155
156
157
158
159
160
            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])

161
162
            mflat.header[CONFIG['DRS_VERSION'][0]] = (CONFIG['VERSION'], CONFIG['DRS_VERSION'][1])

Andrea Bignamini's avatar
Andrea Bignamini committed
163
164
165
            # Add metadata to header
            mflat.header = metadata.add_metadata(mflat.header)

Monica Rainer's avatar
Monica Rainer committed
166
167
168
169
170
171
172
            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))))

173
174
175
176
177
178
            shiftY = varie.shiftY(mflat.data)

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

            mflat = None

Monica Rainer's avatar
Monica Rainer committed
179
180
181
182
183
        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')
184
185
186
187
188
189
190
191
192
193
194
195
196
197
                    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
198
199
200
                except:
                    self.messages.append('There are no flats in the calibration database.')
                    return False
201
202
203
204
205
206
207

            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)

208
            # Since masterflat has been copied from calibration database
209
            # also str and nor will be copied
210
211
212
            self.copy_straighten()
            return False

Monica Rainer's avatar
Monica Rainer committed
213
214
215
        return True


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
    def copy_straighten(self):
        """
        The masterflat has been copied from calibration database
        This function also copy str and nor file from calibration database
        """

        # Get masterflat name from calibdb
        # flat is a fullpath name
        flat = db.extract_dbfile(self.dbconn,'flat')

        # Generate str and nor fullpath name from calibdb
        straight = flat.replace('.fits','_str.fits')
        nor = flat.replace('.fits','_nor.fits')

        # Insert str and nor in calibdb for this night
        db.insert_dbfile(self.dbconn,'flatstr',straight)
        db.insert_dbfile(self.dbconn,'flatnor',nor)

        # Copy str and nor to RED_CALIB directory of current night
        red_straight = os.path.join(CONFIG['RED_CALIB'],os.path.basename(straight))
        shutil.copyfile(straight,red_straight)
        red_nor = os.path.join(CONFIG['RED_CALIB'],os.path.basename(nor))
        shutil.copyfile(nor,red_nor)

        return


Monica Rainer's avatar
Monica Rainer committed
243
244
245
246
247
248
249
250
251
252
253
    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
254
        #extracted = flat.replace('.fits','_ext.fits')
Monica Rainer's avatar
Monica Rainer committed
255

256
257
258
259
260
261
262
263
264
265
266
267
        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
268

269
270
271
272
273
274
275
        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)
276

277
278
279
            shiftY = [''.join(('DY=',str(shift - CONFIG['SHIFT_Y'])))]
            #print shiftY
            args.extend(shiftY)
280

281
        #print args
Monica Rainer's avatar
Monica Rainer committed
282
283
284
285
286
287
        subprocess.call(args)

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

        mflat = ccdproc.CCDData.read(straight, unit=u.adu)
Andrea Bignamini's avatar
Andrea Bignamini committed
288
289
290
291
292
293
294
295
296
297
298
299

        # Update FILENAME in header then
        # add metadata to header and save straight file
        mflat.header[CONFIG['KEYS']['FILENAME']] = os.path.basename(straight)
        mflat.header = metadata.add_metadata(mflat.header)

        hdu = fits.PrimaryHDU(data=mflat.data, header=mflat.header)
        str_fits = fits.HDUList([hdu])
        str_fits.writeto(straight, overwrite=True)



Monica Rainer's avatar
Monica Rainer committed
300
301
302
303
304
305
306
307
308
309
        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
310
311
        #ordini = np.zeros((CONFIG['N_ORD'],CONFIG['XCCD']))

Monica Rainer's avatar
Monica Rainer committed
312
        for x in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
313
314
315
316

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

Monica Rainer's avatar
Monica Rainer committed
317
            # order limits
Monica Rainer's avatar
Monica Rainer committed
318
319
            start = x*CONFIG['W_ORD']
            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
Monica Rainer's avatar
Monica Rainer committed
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338

            # 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)
339
340
                    #norflat[rrow] = fdata[rrow]/pfit
                    norflat[rrow] = np.true_divide(fdata[rrow],pfit)
Monica Rainer's avatar
Monica Rainer committed
341
342
343
344
345
346
                    maskC[rrow] = 0
                    #if not i:
                    #    i = rrow
                    #width = width + 1


Monica Rainer's avatar
Monica Rainer committed
347
348
349
350
351
352
353
354
355
356
        # 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
357
358
359
        fnor = ccdproc.CCDData(norflat, unit=u.adu)
        fnor.header = mflat.header
        #fn = fnor.to_hdu()
Andrea Bignamini's avatar
Andrea Bignamini committed
360
361
362
363
364
365
366
367

        # Update FILENAME in header then
        # add metadata to header and save nor file
        fnor.header[CONFIG['KEYS']['FILENAME']] = os.path.basename(nor)
        fnor.header = metadata.add_metadata(fnor.header)


        hdu = fits.PrimaryHDU(data=fnor.data, header=fnor.header)
Monica Rainer's avatar
Monica Rainer committed
368
        fn = fits.HDUList([hdu])
Andrea Bignamini's avatar
Andrea Bignamini committed
369
        fn.writeto(nor, overwrite=True)
Monica Rainer's avatar
Monica Rainer committed
370

Andrea Bignamini's avatar
Andrea Bignamini committed
371
372
373
        # Copy straight and nor to RED_CALIB directory
        red_straight = os.path.join(CONFIG['RED_CALIB'],os.path.basename(straight))
        shutil.copyfile(straight,red_straight)
Monica Rainer's avatar
Monica Rainer committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
        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