Coding is fun

It’s been almost three weeks since my last post, so I’ve got some catching up to do. The key aspects will follow in the next few posts, but first here is a quick self-promoting “yay!” for myself and my improving Python scripting abilities. In my last week in the office before heading on an assortment of adventures (see following posts), I managed to code up some useful Python modules for use on the pt5m. Hurray! Now I actually deserve at least a tiny bit of credit for being part of the team. Anyway, here is some code to automatically create bias and flat field frames (with assistance from Stu Littlefair). Enjoy:

#!/usr/bin/env python
from pyraf import iraf
import sys
import os, glob
from optparse import OptionParser
import datetime

def getOktal():
    # import current conditions into dictionary 
    # (function parseConditions from home/slodar/scripts/qsi_daemon.py)   
    f = open('/home/slodar/config/conditions_monitor')
    line = f.read().strip()
    condDict = {}
    columns = line.split()
    condDict['cond-date'] = columns[1]+"T"+columns[2]
    dt = datetime.datetime.strptime(condDict['cond-date'],'%Y-%m-%dT%H:%M:%S')
    age = datetime.datetime.utcnow()-dt
    age_seconds = 86400*age.days + age.seconds
    if (age_seconds > 900):
        raise Exception('Conditions log out of date')
    condDict.update( dict([columns[i+3].split('=') for i in range(len(columns)-3)]) )
    # read new paramter skyT as float
    skyT = float(condDict['sky_temp'])
    # convert to oktal 
    # (function temp2oktal from home/slodar/scripts/qsi_daemon.py)   
    p0 = -16.749332
    p1 = 2.88788
    oktal = (skyT-p0)/p1
    return oktal

def cleanupFiles(data_dir):
    cals_dir = data_dir + 'calibs/'
    files=glob.glob(data_dir+'db*')
    files.extend( glob.glob('*.lis') )
    files.extend( glob.glob(cals_dir + '*flat.fits') )
    for filename in files:
        os.remove(filename)
    os.remove('allFits')
    return

def parseCmdLineArgs():
    usage = "usage: %prog [--live] dir_name"
    parser = OptionParser(usage=usage)
    parser.add_option("-l", "--live",action="store_true", 
        dest="live_flag",default=False,help="script being run live on ron?")

    (options, args) = parser.parse_args()  

    live_flag = options.live_flag

    if(len(args) == 1):
        data_dir = args[0]
        # check if data_dir ends in '/', if not, add it
        if not data_dir.endswith('/'): 
            data_dir=data_dir+'/'       
    else:
        # if data directory left blank, use this directory
        data_dir = './'
    return (data_dir,live_flag)

def makeMasterFlatFrames(cals_dir,filt):

    #assume flats are 1x1 (they should be), check for suitable bias
    if not os.path.exists(cals_dir+'bias1x1.fits'):
        print 'no 1x1 bias available this night, aborting'
        return
    else:
        biasFrame = cals_dir+'bias1x1.fits'
    
    matchExpr = "TYPE='SKY' && FILTER='%s'" % filt
    flatlis = iraf.hsel(images='@allFits',fields='$I',expr=matchExpr,Stdout=1)

    #write list files for each filter
    listFile = filt+'flat.lis'
    f = open(listFile, 'w')
    for flat in flatlis:
        f.write(flat + '\n')
    f.close()

    #check if >2 flats exist, if not, continue with next filter    
    if len(flatlis) >2:

        #write de-biassed list filenames for each filter
        dbListFile = filt+'flatdb.lis'
        f = open(dbListFile, 'w')
        for flat in flatlis:    
            head,tail=os.path.split(flat)
            f.write(head+'/db'+tail+'\n')
        f.close()

        #subtract 1x1 bias frame from flats
        iraf.imarith(operand1='@'+listFile,op='-',operand2=biasFrame,result='@'+dbListFile)
        print filt+'-band flats de-biassed'

        # combine de-biassed flat fields, with scaling
        outIm = cals_dir+filt+'flat.fits'
        iraf.imcombine(input='@'+dbListFile,output=outIm,combine='median',reject='avsigclip',scale='mode',statsec='[900:1000,650:750]',Stdout=1)
        print 'Master '+filt+'-band flat created'

        #calculate normalisation constant
        norm_const=iraf.imstat(image=outIm,fields='mean',format=0,Stdout=1)
        norm_const=float(norm_const[0])

        #normalise!        
        iraf.imarith(operand1=outIm,op='/',operand2=norm_const,result=cals_dir+filt+'flatnorm1x1.fits')        
        print 'Master '+filt+'-band flat normalised'

        #create binned flats for each filter     ***using 'average' rather than sum***
        iraf.blkavg(input=cals_dir+filt+'flatnorm1x1.fits',output=cals_dir+filt+'flatnorm2x2.fits',b1=2,b2=2,option='average')
        print 'Master '+filt+'-band 2x2 flat created'

        iraf.blkavg(input=cals_dir+filt+'flatnorm1x1.fits',output=cals_dir+filt+'flatnorm3x3.fits',b1=3,b2=3,option='average')
        print 'Master '+filt+'-band 3x3 flat created'

    # if no flats of this filter exist, move to next filter
    else:
        print 'Less than 3 '+filt+'-band flats available. Filter skipped'
        return


def makeMasterBiasFrames(cals_dir,bin):
    bin = str(bin)
    
    matchExpr = "TYPE='BIAS' && XBINNING=%s" % bin
    
    # find matching biasses
    biaslis = iraf.hsel(images='@allFits',fields='$I',expr=matchExpr,Stdout=1)

    #write list files for each 
    listFile = 'bias'+bin+'.lis'
    f = open(listFile, 'w')
    for bias in biaslis:
        f.write(bias + '\n')
    f.close()
    
    #combine to create master biases
    if len(biaslis) > 2:
        imOut = cals_dir + 'bias'+bin+'x'+bin+'.fits'
        iraf.imcombine(input='@'+listFile,output=imOut,combine='median',Stdout=1)
        print 'Binning '+bin+'x'+bin+' bias created'
    else:
        print 'Too few '+bin+'x'+bin+' biases to combine'


if __name__ == "__main__":
    data_dir, live_flag = parseCmdLineArgs()

    # allow iraf to overwrite existing files
    iraf.reset(clobber='yes')

    # always put calibrations in calibs, for consistency
    cals_dir=data_dir+'calibs/'
    if not os.path.exists(cals_dir):
        os.mkdir(cals_dir)

    # 'hsel' is not always available by default, load images package
    iraf.images(_doprint=0)

    # select all frames
    # filename template must ignore any _cat.fits files...
    iraf.files(data_dir+'r*[0-9].fits',Stdout='allFits')


    ###############################
    ## Biasses ##
    for bin in ['1','2','3']:
        try:
            makeMasterBiasFrames(cals_dir,bin)
        except Exception as e:
            print 'Encountered error making bias'
            print e
            print 'Continuing'
            

    # when running live, check for cloud levels
    # abort script at this point if cloud present
    if live_flag:
        try:
            oktal = getOktal()
        except Exception as e:
            print 'could not get cloud conditions, aborting'
            print e
            cleanupFiles(data_dir)
            sys.exit()

        # if oktal above zero, abort script here
        if oktal > 0.0:
            print 'Suspect cloudy conditions. Aborting script here. No master flats created.'
            ## delete used files ##
            cleanupFiles(data_dir)
            sys.exit()


    ################################
    ## Flats ##
    filters = ['H','B','V','R','I']
    for filt in filters:
        try:
            makeMasterFlatFrames(cals_dir,filt)
        except Exception as e:
            print 'Encountered error making flat'
            print e
            print 'Continuing'
        
    cleanupFiles(data_dir)
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s