flatframes.py 9.81 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
54
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
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
135
"""
Last modified: 2017-03-07

Reduction of the flat-fields:
- check the exposure time
- 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
            else:
                self.flatlist.append(flat)
        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)

            combineflat = ccdproc.Combiner(self.flatcorr)

            # mask the pixels using ron/gain
Monica Rainer's avatar
Monica Rainer committed
136
137
            #combineflat.sigma_clipping(func=np.ma.mean, dev_func=varie.stdcombine)
            combineflat.sigma_clipping(func=np.ma.median, dev_func=varie.stdcombine)
Monica Rainer's avatar
Monica Rainer committed
138
            mflat = combineflat.average_combine()
Monica Rainer's avatar
Monica Rainer committed
139

Monica Rainer's avatar
Monica Rainer committed
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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
            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
        MODIFIED: no normalization, only straightened and mask C creation
        """
        # straighten the flat

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

        straight = flat.replace('.fits','_str.fits')
        nor = flat.replace('.fits','_nor.fits')

        args = [CONFIG['STRAIGHT'],flat,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('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
206
        for x in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
207
            # order limits
Monica Rainer's avatar
Monica Rainer committed
208
209
            start = x*CONFIG['W_ORD']
            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
Monica Rainer's avatar
Monica Rainer committed
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

            # 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


        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