casa: task_sdbaseline.py Source File

Go to the documentation of this file.
00001 import os
00002 from taskinit import *
00003 
00004 import sdutil
00005 import asap as sd
00006 from asap._asap import Scantable
00007 from asap import _to_list
00008 from asap.scantable import is_scantable
00009 import math
00010 import pylab as pl
00011 
00012 def sdbaseline(infile, antenna, fluxunit, telescopeparm, specunit, restfreq, frame, doppler, scanlist, field, iflist, pollist, tau, masklist, maskmode, thresh, avg_limit, edge, blfunc, order, npiece, applyfft, fftmethod, fftthresh, addwn, rejwn, clipthresh, clipniter, verify, verbose, bloutput, blformat, showprogress, minnrow, outfile, outform, overwrite, plotlevel):
00013     
00014     casalog.origin('sdbaseline')
00015     
00016     try:
00017         worker = sdbaseline_worker(**locals())
00018         worker.initialize()
00019         worker.execute()
00020         worker.finalize()
00021         
00022     except Exception, instance:
00023         sdutil.process_exception(instance)
00024         raise Exception, instance
00025 
00026 class sdbaseline_worker(sdutil.sdtask_template):
00027     def __init__(self, **kwargs):
00028         super(sdbaseline_worker,self).__init__(**kwargs)
00029         self.suffix = '_bs'
00030 
00031     def parameter_check(self):
00032         if (self.order < 0) and (self.blfunc in ['poly', 'chevishev']):
00033             msg = 'Negative order of baseline polynomial given. Exit without baselining.'
00034             casalog.post(msg, priority = 'SEVERE')
00035             raise Exception(msg)
00036 
00037         if (self.npiece <= 0) and (self.blfunc == 'cspline'):
00038             msg = 'The parameter npiece must be greater than 0. Exit without baselining.'
00039             casalog.post(msg, priority='SEVERE')
00040             raise Exception(msg)            
00041 
00042     def initialize_scan(self):
00043         sorg = sd.scantable(self.infile, average=False, antenna=self.antenna)
00044         
00045         sel = self.get_selector()
00046         sorg.set_selection(sel)
00047         del sel
00048 
00049         # Copy scantable when usign disk storage not to modify
00050         # the original table.
00051         if is_scantable(self.infile) and self.is_disk_storage:
00052             self.scan = sorg.copy()
00053         else:
00054             self.scan = sorg
00055         del sorg
00056 
00057     def execute(self):
00058         engine = sdbaseline_engine(self)
00059         engine.prologue()
00060 
00061         # check if the data contains spectra
00062         if (self.scan.nchan()==1):
00063            s = "Cannot process the input data. It contains only single channel data."
00064            raise Exception, s
00065    
00066         # set various attributes to self.scan
00067         self.set_to_scan()
00068         
00069         scanns = self.scan.getscannos()
00070         sn=list(scanns)
00071         casalog.post( "Number of scans to be processed: %d" % (len(sn)) )
00072         
00073         # Warning for multi-IF data
00074         #if ( len(s.getifnos()) > 1 and not maskmode == 'auto' ):
00075         if ( len(self.scan.getifnos()) > 1 and isinstance(self.masklist,list) and not self.maskmode == 'auto' ):
00076             casalog.post( 'The scantable contains multiple IF data.', priority = 'WARN' )
00077             casalog.post( 'Note the same mask(s) are applied to all IFs based on CHANNELS.', priority = 'WARN' )
00078             casalog.post( 'Baseline ranges may be incorrect for all but IF=%d.' % (self.scan.getif(0)), priority = 'WARN' )
00079             
00080         # do opacity (atmospheric optical depth) correction
00081         sdutil.doopacity(self.scan, self.tau)
00082     
00083         if (self.order < 0):
00084             casalog.post('Negative order of baseline polynomial given. Exit without baselining.', priority = 'WARN')
00085             return
00086         engine.drive()
00087 
00088         # Plot final spectrum
00089         engine.epilogue()
00090 
00091     def save(self):
00092         sdutil.save(self.scan, self.project, self.outform, self.overwrite)
00093 
00094 class sdbaseline_engine(sdutil.sdtask_engine):
00095     def __init__(self, worker):
00096         super(sdbaseline_engine,self).__init__(worker)
00097 
00098         clip_keys = ['clipthresh', 'clipniter']
00099         self.poly_keys = ['order']
00100         self.chevishev_keys = self.poly_keys + clip_keys
00101         self.cspline_keys = ['npiece'] + clip_keys
00102         self.sinusoid_keys = ['applyfft', 'fftmethod', 'fftthresh', 'addwn', 'rejwn'] + clip_keys
00103         self.auto_keys = ['thresh', 'avg_limit', 'edge']
00104         self.list_keys = ['masklist']
00105         self.interact_keys = []
00106 
00107         self.baseline_param = sdutil.parameter_registration(self)
00108         self.baseline_func = None
00109 
00110     def __del__(self):
00111         del self.baseline_param
00112 
00113     def __register(self, key, attr=None):
00114         self.baseline_param.register(key, attr)
00115 
00116     def __get_param(self):
00117         return self.baseline_param.get_registered()
00118 
00119     def prologue(self):
00120         if ( abs(self.plotlevel) > 1 ):
00121             casalog.post( "Initial Raw Scantable:" )
00122             self.worker.scan._summary()
00123 
00124     def drive(self):
00125         scan = self.worker.scan
00126         self.__init_blfile()
00127 
00128         nrow = scan.nrow()
00129 
00130         # parse string masklist
00131         if isinstance(self.masklist,list):
00132                 maskdict = {'': self.masklist}
00133         else:
00134                 maskdict = scan.parse_maskexpr(self.masklist)
00135         basesel = scan.get_selection()
00136 
00137         # configure baseline function and its parameters
00138         self.__configure_baseline()
00139         
00140         for sif, lmask in maskdict.iteritems():
00141             if len(sif) > 0:
00142                 #s.set_selection(selection=(basesel+sd.selector(ifs=[int(sif)])))
00143                 sel = sd.selector(basesel)
00144                 sel.set_ifs([int(sif)])
00145                 scan.set_selection(sel)
00146                 del sel
00147                 msg = "Working on IF%s" % (sif)
00148                 casalog.post(msg)
00149                 if (self.maskmode == 'interact'): print "===%s===" % (msg)
00150                 del msg
00151 
00152             msk = None
00153 
00154             if (self.maskmode == 'interact'):
00155                 msk = sdutil.interactive_mask(scan, lmask, False, purpose='to baseline spectra')
00156                 msks = scan.get_masklist(msk)
00157                 if len(msks) < 1:
00158                     msg = 'No channel is selected. Exit without baselining.'
00159                     casalog.post(msg, priorigy='SEVERE')
00160                     raise Exception(msg)
00161 
00162                 casalog.post( 'final mask list ('+scan._getabcissalabel()+') ='+str(msks) )
00163                 #header += "   Fit Range: "+str(msks)+"\n"
00164                 del msks
00165             else:
00166                 # Use baseline mask for regions to INCLUDE in baseline fit
00167                 # Create mask using list, e.g. masklist=[[500,3500],[5000,7500]]
00168                 if (len(lmask) > 0): msk = scan.create_mask(lmask)
00169 
00170             # register IF dependent mask
00171             self.__register('mask', msk)
00172 
00173             # call target baseline function with appropriate parameter set 
00174             self.baseline_func(**self.__get_param())
00175 
00176 ##             if (maskmode == 'auto'):
00177 ##                 if (blfunc == 'poly'):
00178 ##                     scan.auto_poly_baseline(mask=msk,order=order,edge=edge,threshold=thresh,chan_avg_limit=avg_limit,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00179 ##                 elif (blfunc == 'chebyshev'):
00180 ##                     scan.auto_chebyshev_baseline(mask=msk,order=order,clipthresh=clipthresh,clipniter=clipniter,edge=edge,threshold=thresh,chan_avg_limit=avg_limit,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00181 ##                 elif (blfunc == 'cspline'):
00182 ##                     scan.auto_cspline_baseline(mask=msk,npiece=npiece,clipthresh=clipthresh,clipniter=clipniter,edge=edge,threshold=thresh,chan_avg_limit=avg_limit,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00183 ##                 elif (blfunc == 'sinusoid'):
00184 ##                     scan.auto_sinusoid_baseline(mask=msk,applyfft=applyfft,fftmethod=fftmethod,fftthresh=fftthresh,addwn=addwn,rejwn=rejwn,clipthresh=clipthresh,clipniter=clipniter,edge=edge,threshold=thresh,chan_avg_limit=avg_limit,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00185 ##             else:
00186 ##                 if (blfunc == 'poly'):
00187 ##                     scan.poly_baseline(mask=msk,order=order,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00188 ##                 elif (blfunc == 'chebyshev'):
00189 ##                     scan.chebyshev_baseline(mask=msk,order=order,clipthresh=clipthresh,clipniter=clipniter,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00190 ##                 elif (blfunc == 'cspline'):
00191 ##                     scan.cspline_baseline(mask=msk,npiece=npiece,clipthresh=clipthresh,clipniter=clipniter,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00192 ##                 elif (blfunc == 'sinusoid'):
00193 ##                     scan.sinusoid_baseline(mask=msk,applyfft=applyfft,fftmethod=fftmethod,fftthresh=fftthresh,addwn=addwn,rejwn=rejwn,clipthresh=clipthresh,clipniter=clipniter,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile,csvformat=csvformat,insitu=True)
00194 
00195                 # the above 14 lines will eventually shrink into the following 2 commands:
00196                 #
00197                 #sbinfo = s.create_sbinfo(blfunc=blfunc,order=order,npiece=npiece,applyfft=applyfft,fftmethod=fftmethod,fftthresh=fftthresh,addwn=addwn,rejwn=rejwn,\
00198                 #                         masklist=masklist,maskmode=maskmode,edge=edge,threshold=threshold,chan_avg_limit=chan_avg_limit,\
00199                 #                         clipthresh=clipthresh,clipniter=clipniter)
00200                 #s.sub_baseline(sbinfo=sbinfo,plot=verify,showprogress=showprogress,minnrow=minnrow,outlog=verbose,blfile=blfile)
00201                 #
00202                 # where
00203                 # sbinfo = {'func':funcinfo, 'mask':maskinfo, 'clip':clipinfo}
00204                 # and
00205                 # funcinfo should be one of the follows:
00206                 #     funcinfo = {'type':'poly', 'params':{'order':order}}
00207                 #     funcinfo = {'type':'cspline', 'params':{'npiece':npiece}}
00208                 #     funcinfo = {'type':'sspline', 'params':{'lambda':lambda}}
00209                 #     funcinfo = {'type':'sinusoid', 'params':{'applyfft':applyfft, 'fftmethod':fftmethod, 'fftthresh':fftthresh, 'addwn':addwn, 'rejwn':rejwn}}
00210                 # maskinfo should be one of the follows:
00211                 #     maskinfo = {'base':masklist, 'aux':{'type':'auto', 'params':{'edge':edge, 'threshold':thresh, 'chan_avg_limit':avg_limit}}}
00212                 #     maskinfo = {'base':masklist, 'aux':{'type':'list'}}
00213                 #     maskinfo = {'base':masklist, 'aux':{'type':'interactive'}}
00214                 # clipinfo should be:
00215                 #     clipinfo = {'clipthresh':clipthresh, 'clipniter':clipniter}
00216 
00217             # reset selection
00218             if len(sif) > 0: scan.set_selection(basesel)
00219             
00220     def epilogue(self):
00221         if ( abs(self.plotlevel) > 0 ):
00222             pltfile=self.project+'_bsspec.eps'
00223             sdutil.plot_scantable(self.worker.scan, pltfile, self.plotlevel)
00224 
00225     def __configure_baseline(self):
00226         # determine what baseline function should be called
00227         funcname = '%s_baseline'%(getattr(self,'blfunc').lower())
00228         if self.maskmode.lower() == 'auto':
00229             funcname = 'auto_' + funcname
00230         self.baseline_func = getattr(self.worker.scan,funcname)
00231 
00232         # register parameters for baseline function
00233         # parameters for auto mode
00234         if self.maskmode.lower() == 'auto':
00235             self.__register('threshold', 'thresh')
00236             self.__register('chan_avg_limit', 'avg_limit')
00237             self.__register('edge')
00238 
00239         # parameters that depends on baseline function
00240         keys = getattr(self, '%s_keys'%(self.blfunc.lower()))
00241         for k in keys:
00242             self.__register(k)
00243 
00244         # common parameters
00245         self.__register('plot', 'verify')
00246         self.__register('showprogress')
00247         self.__register('minnrow')
00248         self.__register('outlog', 'verbose')
00249         self.__register('blfile')
00250         self.__register('csvformat', (self.blformat.lower() == "csv"))
00251         self.__register('insitu', True)
00252 
00253 ##         for (k,v) in self.__get_param().items():
00254 ##             casalog.post('%s=%s'%(k,v),'WARN')
00255 
00256     def __init_blfile(self):
00257         if self.bloutput:
00258             self.blfile = self.project + "_blparam.txt"
00259 
00260             if (self.blformat.lower() != "csv"):
00261                 f = open(self.blfile, "w")
00262 
00263                 blf = self.blfunc.lower()
00264                 mm = self.maskmode.lower()
00265                 if blf == 'poly':
00266                     ftitles = ['Fit order']
00267                 elif blf == 'chevishev':
00268                     ftitles = ['Fit order', 'clipThresh', 'clipNIter']
00269                 elif blf == 'cspline':
00270                     ftitles = ['nPiece', 'clipThresh', 'clipNIter']
00271                 else: # sinusoid
00272                     ftitles = ['applyFFT', 'fftMethod', 'fftThresh', 'addWaveN', 'rejWaveN', 'clipThresh', 'clipNIter']
00273                 if mm == 'auto':
00274                     mtitles = ['Threshold', 'avg_limit', 'Edge']
00275                 elif mm == 'list':
00276                     mtitles = ['Fit Range']
00277                 else: # interact
00278                     mtitles = []
00279                     
00280                 fkeys = getattr(self, '%s_keys'%(self.blfunc))
00281                 mkeys = getattr(self, '%s_keys'%(self.maskmode))
00282 
00283                 info = [['Source Table', self.infile],
00284                         ['Output File', self.project],
00285                         ['Flux Unit', self.worker.scan.get_fluxunit()],
00286                         ['Abscissa', self.worker.scan.get_unit()],
00287                         ['Function', self.blfunc]]
00288                 for i in xrange(len(ftitles)):
00289                     info.append([ftitles[i],getattr(self,fkeys[i])])
00290                 info.append(['Mask mode', self.maskmode])
00291                 for i in xrange(len(mtitles)):
00292                     info.append([mtitles[i],getattr(self,mkeys[i])])
00293 
00294                 separator = "#"*60 + "\n"
00295                 
00296                 f.write(separator)
00297                 for i in xrange(len(info)):
00298                     f.write('%12s: %s\n'%tuple(info[i]))
00299                 f.write(separator)
00300                 f.close()
00301         else:
00302             self.blfile = ""        
00303