from darepype.drp import DataFits # pipeline data object class from darepype.drp.stepmiparent import StepMIParent # pipestep Multi-Input parent from darepype.tools.steploadaux import StepLoadAux # pipestep steploadaux object class from astropy.io import fits #package to recognize FITS files import numpy as np import logging from astropy.stats import mad_std # The median absolute deviation, a more robust estimator than std. from astropy.time import Time #?import sys import os #?import astropy #? from astropy import units as u class StepMasterFlatCCD(StepLoadAux, StepMIParent): ''' Pipe step object to create master flats for HDR images taken with SBIG CMOS camera. Inputs are sets of high-gain (bin1H), and low-gain (bin1L) RAW flat images, all taken through the same spectral filter and electronic gain. Required auxiliary files include two PFIT files, one each for high and low gain data. PFIT files contain data required to synthesize dark images exactly matching the exposure times of the input flat images. PFIT files contain a 3D image with a dark current image in the first plane and a bias image in the second plane. Output is a single file that contains a two-plane 3D imgage in the primary HDU. The first plane is a high-gain flat and the second plane is a low-gain flat. A second HDU contains a relative gain image (the ratio of the electronic gains of the pixels in the high and low-gain images. A third HDU contains a table comprising statistical information about the input images and instrumental and environmental conditions during the exposures. ''' def __init__(self): """ Constructor: Initialize data objects and variables """ # Call superclass constructor (calls setup) super(StepMasterFlatHdr,self).__init__() self.flatloaded = False # Indicates if flat file has been loaded self.flat = None # Numpy array object containing flat values self.flatname = '' # name of selected flat file self.biasloaded = False # Indicates if flat file has been loaded self.bias = None # Numpy array object containing flat values self.biasname = '' # name of selected flat file self.darkloaded = False # Indicates if flat file has been loaded self.dark = None # Numpy array object containing flat values self.darkname = '' # name of selected flat file # Finish up. self.log.debug('Init: done') def setup(self): ''' Names and Parameters need to be set here Sets the internal names for the function and for saved files. Defines the input parameters for the current pipe step. Setup() is called at the end of __init__ The parameters are stored in a list of entries containing the following information: - name: The name for the parameter. This name is used when calling the pipe step from command line or python shell. It is also used to identify the parameter in the pipeline configuration file. - default: A default value for the parameter. If nothing, set '' for strings, 0 for integers and 0.0 for floats - help: A short description of the parameter. ''' ### Set Names # Name of the pipeline reduction step self.name='masterflathdr' # Shortcut for pipeline reduction step and identifier for saved file names. self.procname = 'HDRFLAT' ### Set Logger for this pipe step self.log = logging.getLogger('pipe.step.%s' % self.name) ### Set Parameter list # Clear Parameter list self.paramlist = [] # Append parameters self.paramlist.append(['combinemethod','median', 'Specifies how the files should be combined - options are median, average, sum']) self.paramlist.append(['outputfolder','', 'Output directory location - default is folder with input files']) self.paramlist.append(['hotpxlim', 99.5, 'Hot pixel limit percentile']) self.paramlist.append(['reload', False,'Set to True to look for new pfit files for every input']) self.paramlist.append(['outputfolder2', '','Alternate output directory path']) self.paramlist.append(['gainpcntlim', 0.3,'gain quality threshold']) self.paramlist.append(['dstdpcntlim', 10.0,'dstd quality threshold']) self.paramlist.append(['numfilelim', 8,'Minumum number of input files']) self.paramlist.append(['print_switch', False,'Set True to turn on print statements']) ### Set parameters for StepLoadAux self.loadauxsetup('MDARK') self.loadauxsetup('MBIAS') def timesortHDR(self, datalist, date_key = 'date-obs'): ''' Sorts a list of fits files by a header keyword with a date/time value. In the case of an SBIG CMOS camera RAW file, the date/time is read from the second HDU. Arguments: datalist = a list of DataFits objects date_key = the header keyword containing the time/date data Returns: tfiles = the sorted file list utime = a list of the unix times of the observations Author(s): Al Harper, Finian Ashmead Modified: 210807, 210815, 220802 (modified for use in pipe step) Version: 1.1 (pipe) ''' date_obs = [] # Make a list to hold the date-obs keyword strings. for d in datalist: if '_bin1L' in d.filename: head = d.getheader(d.imgnames[1]) # Get the header of the second HDU (index = [1]). date_obs.append(head[date_key]) # Add date information to list. of string objects. else: head = d.getheader() # Get the header of the primary HDU (index = [0]). date_obs.append(head[date_key]) # Add date information to list. of string objects. t = Time(date_obs, format='isot', scale='utc') # Make an astropy time object in 'isot' format. tsort = np.argsort(t) # Make a list of indices that will sort by date_obs. tfiles = [] utime = [] for i in tsort: tfiles.append(datalist[i]) utime.append(t[i].unix) return tfiles, utime def run(self): ''' Process sets of high-gain and low-gain RAW flatfiles (datain). This code is based on Al's Jupyter notebook make_flat_HDR_auto_54. ''' ''' Define boolean to enable print statements for debugging. Set = True to print. ''' pt = self.getarg('print_switch') ''' Make two separate lists of input DataFits objects, one for high gain images and one for low gain images. Then sort files by exposure time in order to get matched lists of high and low gain images associated with the same identical exposure. ''' flatfiles = [f for f in os.listdir(flatpath) if '.fit' in f and binning in f \ and filekey in f and specfltr in f and 'RAW.fit' in f and '._' not in f] flatfiles, utimeH = timesortHDR(flatfiles, flatpath, date_key = 'date-obs', \ print_list = True) df = read_oneDF(flatfiles, 0, flatpath) rows, cols = df.image.shape[0], df.image.shape[1] ''' Make a 3D image of flat files. Dimensions are [stack, row, col]. ''' numflats = len(flatfiles) flatheadlist, flatstats = [[], []], [[]] flatimage = np.zeros((numflats, rows, cols)) flatbaseheader = flatheadlist[0] ''' Compile a list of exposure times. Set keyword variables for max and min exposure times for inclusion in header. ''' flatexptimes = [] imedian = flatstats['median'] for j in range(len(flatfiles)): exptime = flatheadlist[j].header['exptime'] flatexptimes.append(exptime) xtimemin = np.min(flatexptimes) xtimemax = np.max(flatexptimes) if pt: print('xtimemin =', xtimemin, ' xtimemax =', xtimemax) ''' Construct a stack of dark images to match the flat image exposure times. ''' darkimage = np.zeros_like(flatimage) for j in range(flatimage.shape[0]): darkimage[j] = ((dark - bias) * (flatexptimes[j])/darkexptime) + bias ''' Subtract interpolated darks from the flat images. ''' flatimageDS = flatimage - darkimage ''' Normalize each of the images in the flat image stack to its own median (after applying hotpix and gain masks). From this point on, there will be nans for every pixel in the hotpix and gain masks. ''' flatimageDSN = np.zeros_like(flatimageDS) flatmediansDS = np.zeros((numflats)) for j in range(numflats): flatimageDS[j][hotpix] = np.nan flatmediansDS[j] = np.nanmedian(flatimageDS[j]) flatimageDSN[j] = flatimageDS[j] / flatmediansDS[j] ''' Make median flats from the stacks and compute their respective median, mean, std, and mad. The result is a 2D image. ''' flat = np.nanmedian(flatimageDSN, axis=0) flatmedian, flatmean, flatstd, flatmadstd = [0], [0], [0], [0] flatmedian = np.nanmedian(flat) # Median of the median flat. flatmean = np.nanmean(flat) # Mean of the median flat. flatstd = np.nanstd(flat) # mad_std of the median flat. flatmadstd = mad_std(flat,ignore_nan=True) # mad_std of the median flat. if pt: print('Median flat median =', flatmedian ) if pt: print('Mean flat median =', flatmean ) if pt: print('Median flat std =', flatstd) if pt: print('Median flat mad_std =', flatmadstd) if pt: print('Shape of flatimgage:', flatimage.shape) if pt: print('Shape of flat:', flat.shape) if pt: print('') ''' if pt: print a list of the medians of the dark-subtracted high and low gain flat images. ''' if pt: print('List of medians of dark-subtracted flat images') for i in range(numflats): if pt: print('{:<3}{:<10.2f}'.format(i, flatmediansDS[i])) if pt: print('') ''' Find minimum and maximum values of the image medians for inclusion in output header. ''' medmin, medmax = [0], [0] medmin = np.nanmin(flatmediansDS) medmax = np.nanmax(flatmediansDS) if pt: print('minimum medians =', medmin) if pt: print('maximum medians =', medmax) if pt: print('') ''' Re-normalize the median flats (H and L) to their medians and mask out pixels with gains greater than or less than 1.0 from the median gain. Replace the masked pixels with np.nan. Hence, when a sky image is divided by the flat, those pixels will also be masked (will be np.nan) in the flat-fielded sky image (in addition to any pixela aleady replaced with nans using the hotpix mask). ''' mflat = np.zeros_like(flat) mflat = flat / flatmedian mflatmedian = np.nanmedian(mflat) mflatmadstd = mad_std(mflat, ignore_nan=True) if pt: print('mflatmedian, mflatmadstd =', mflatmedian, mflatmadstd) if pt: print('') ''' Compute the differences of each of the individual normalized flat images from the median flat of the stack. Compute the median, std, and mad_std of each of the difference images. Define some statistical measures of the uniformity of the difference images that can be reported in header keywords. ''' vmx = mflatmadstd * 2.0 vmn = - mflatmadstd * 2.0 difimage = np.zeros_like(flatimageDS) # Ratio of inidvidual flat images to the median flat image. difmadstd = np.zeros((numflats)) difstd = np.zeros((numflats)) difmean = np.zeros((numflats)) difmedian = np.zeros((numflats)) for i in range(flatimageDS.shape[0]): difimage[i] = flatimageDSN[i]-mflat difmadstd[i] = mad_std(difimage[i],ignore_nan=True) difstd[i] = np.nanstd(difimage[i]) difmean[i] = np.nanmean(difimage[i]) difmedian[i] = np.nanmedian(difimage[i]) dstdmean, dstdmax, dstdmin = np.nanmean(difstd), np.nanmax(difstd), np.nanmin(difstd) dstdpcnt = (dstdmax - dstdmin)*100/dstdmean ''' Create output DataFits object and fill with image data. ''' self.dataout = DataFits(config=self.config) self.dataout.header = flatbaseheader.copy() self.dataout.image = mflat ''' Create output file name. ''' infolder = os.path.split(flatfiles.filename) lfname = os.path.split(len(flatfiles)-1).filename) ffname = os.path.split(flatfiles[0].filename) if pt: print('lfname =', lfname) if pt: print('ffname =', ffname) lf = lfname.split('_') # Last filename of time-sorted list. ff = ffname.split('_') # First filename of time-sorted list. flatname = 'mflat_'+ff[1]+'_'+ff[3]+'DR'+'_'+ff[4]+'_'+ff[5]+'-'+lf[5]+'_'+ff[6]+'_'+ff[7]+'_'+'XXX.fits' if pt: print('flatname =',flatname) if pt: print('') ''' Rename output filename ''' gainpcntlim = self.getarg('gainpcntlim') dstdpcntlim = self.getarg('dstdpcntlim') numfilelim = self.getarg('numfilelim') quality = gainpcnt < gainpcntlim and dstdpcnt < dstdpcntlim and numflats >= numfilelim outputfolder = self.getarg('outputfolder') outputfolder2 = self.getarg('outputfolder2') if (outputfolder != '') and (quality == True): outputfolder = os.path.expandvars(outputfolder) self.dataout.filename = os.path.join(outputfolder, flatname) elif (outputfolder != '') and (quality == False): outputfolder = os.path.expandvars(outputfolder2) self.dataout.filename = os.path.join(outputfolder2, flatname) else: self.dataout.filename = os.path.join(infolder, flatname) if pt: print('outputfolder = ',outputfolder) if pt: print('outputfolder2 = ',outputfolder2) if pt: print('dstdpcntlim = ',dstdpcntlim) if pt: print('gainpcntlim = ',gainpcntlim) if pt: print('numfilelim = ',numfilelim) ''' Create numpy arrays with information about input data from headers and timesortDF functions. ''' # From statistical calculations on input data: imedianH, imadH, imeanH, istdH = flatstats['median'], flatstats['mad']\ , flatstats['mean'], flatstats['std'] # From time_sortDF: utime = np.asarray(utimeH) etime = utime - utime[0] # Time elepased from beginning of first exposure of sequence. # From header: ambient = np.zeros((numflats)) primary = np.zeros((numflats)) secondar = np.zeros((numflats)) dewtem1 = np.zeros((numflats)) for i in range(numflats): ambient[i] = flatheadlist[i].header['ambient'] primary[i] = flatheadlist[i].header['primary'] secondar[i] = flatheadlist[i].header['secondar'] dewtem1[i] = flatheadlist[i].header['dewtem1'] # Make a column with the file sequence numbers. index = np.arange(numflats) ''' Now put derived and header data into a fits table and add it to the output object. ''' ''' Make file identifiers. The file identifiers are the date and time fields of the input filenames. ''' IDs = [] for i in range(numflats): fi = highgainlist[i].filename.split('_') IDs.append(fi[4]+'_'+fi[5]) fileIDs = np.asarray(IDs) ''' Make a list of fits column objects. ''' tcols = [] tcols.append(fits.Column(name='index', format='I', array=index)) tcols.append(fits.Column(name='fileID', format='20A', array=fileIDs)) tcols.append(fits.Column(name='median', format='D', array=imedianH, unit='ADU')) tcols.append(fits.Column(name='mean', format='D', array=imeanH, unit='ADU')) tcols.append(fits.Column(name='std', format='D', array=istdH, unit='ADU')) tcols.append(fits.Column(name='mad', format='D', array=imadH, unit='ADU')) tcols.append(fits.Column(name='dmedian', format='D', array=difmedian, unit='ADU')) tcols.append(fits.Column(name='dmean', format='D', array=difmean, unit='ADU')) tcols.append(fits.Column(name='dstd', format='D', array=difstd, unit='ADU')) tcols.append(fits.Column(name='dmad', format='D', array=difmadstd, unit='ADU')) tcols.append(fits.Column(name='ambient', format='D', array=ambient, unit='C')) tcols.append(fits.Column(name='primary', format='D', array=primary, unit='C')) tcols.append(fits.Column(name='secondar', format='D', array=secondar, unit='C')) tcols.append(fits.Column(name='dewtem1', format='D', array=dewtem1, unit='C')) tcols.append(fits.Column(name='elapsed time', format='D', array=etime, unit='seconds')) ''' Make table and add it to the output object. ''' c = fits.ColDefs(tcols) table = fits.BinTableHDU.from_columns(c) tabhead = table.header self.dataout.tableset(table.data, tablename = 'table', tableheader=tabhead) ''' Populate output header with new keyword data. ''' self.dataout.header['notes'] = '1st HDU: 2D mflat image' self.dataout.header['notes3'] = 'Table HDU: statistical and environmental data' self.dataout.header['imagetyp'] = 'MFLAT' self.dataout.header['bzero'] = 0.0 self.dataout.header['ambient'] = np.nanmean(ambient) self.dataout.header['primary'] = np.nanmean(primary) self.dataout.header['secondar'] = np.nanmean(secondar) self.dataout.header['dewtem1'] = np.nanmean(dewtem1) self.dataout.setheadval('xtimemin', xtimemin, 'Minimum exposure in the set of flats') self.dataout.setheadval('xtimemax', xtimemax, 'Maximum exposure in the set of flats') self.dataout.setheadval('medminH', medmin, 'Minimum high gain median in the set of flats') * self.dataout.setheadval('medmaxH', medmax, 'Maximum high gain median in the set of flats') self.dataout.setheadval('dstdmean', dstdmean, 'Mean std in the set of difference images') self.dataout.setheadval('dstdpcnt', dstdpcnt, '(dstdmax-dstdmin)*100.0/dstdmean') self.dataout.setheadval('gainmean', grat_mean, 'Mean gain ratio') self.dataout.setheadval('gainpcnt', gainpcnt, '(gainmax-gainmin)*100.0/gainmean') self.dataout.setheadval('numfiles', numflats, 'Number of RAW exposures in input datasets') self.dataout.setheadval('hotpxlim', hotpxlim, 'Upper limit percentile for unmasked dark current') ''' Add a history keyword ''' self.dataout.setheadval('HISTORY','HDR Master Flat: made from %d x 2 files' % numflats) if pt: print(self.dataout.header) if __name__ == '__main__': """ Main function to run the pipe step from command line on a file. Command: python stepparent.py input.fits -arg1 -arg2 . . . Standard arguments: --config=ConfigFilePathName.txt : name of the configuration file -t, --test : runs the functionality test i.e. pipestep.test() --loglevel=LEVEL : configures the logging output for a particular level -h, --help : Returns a list of """ StepMasterFlatCCD().execute()