flatframes.py 10.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
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
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
        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.
        """
        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:

                # subdtract the masterdak
                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
135
            if len(self.flatcorr) > 1:
Monica Rainer's avatar
Monica Rainer committed
136

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

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

            mflat = None

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

        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')
                except:
                    self.messages.append('There are no flats in the calibration database.')
                    return False
        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
190
        #extracted = flat.replace('.fits','_ext.fits')
Monica Rainer's avatar
Monica Rainer committed
191

192
193
194
        #args = [CONFIG['STRAIGHT'],flat,straight,CONFIG['STRAIGHT_OPT']]
        args = [CONFIG['STRAIGHT'],flat,straight]
        args.extend(CONFIG['STRAIGHT_OPT'])
Monica Rainer's avatar
Monica Rainer committed
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
        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
214
215
        #ordini = np.zeros((CONFIG['N_ORD'],CONFIG['XCCD']))

Monica Rainer's avatar
Monica Rainer committed
216
        for x in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
217
218
219
220

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

Monica Rainer's avatar
Monica Rainer committed
221
            # order limits
Monica Rainer's avatar
Monica Rainer committed
222
223
            start = x*CONFIG['W_ORD']
            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
Monica Rainer's avatar
Monica Rainer committed
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

            # 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
250
251
252
253
254
255
256
257
258
259
        # 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
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
        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