casa: task_sdplot.py Source File

Go to the documentation of this file.
00001 import os
00002 from taskinit import *
00003 
00004 import asap as sd
00005 import pylab as pl
00006 from asap import _to_list
00007 from asap.scantable import is_scantable
00008 import sdutil
00009 
00010 def sdplot(infile, antenna, fluxunit, telescopeparm, specunit, restfreq, frame, doppler, scanlist, field, iflist, pollist, beamlist, scanaverage, timeaverage, tweight, polaverage, pweight, kernel, kwidth, plottype, stack, panel, flrange, sprange, linecat, linedop, subplot, colormap, linestyles, linewidth, histogram, center, cell, header, headsize, plotstyle, margin, legendloc, outfile, overwrite):
00011 
00012     casalog.origin('sdplot')
00013 
00014     ###
00015     ### Now the actual task code
00016     ###
00017     try:
00018         # Plotting
00019         worker = sdplot_worker(**locals())
00020         worker.initialize()
00021         worker.execute()
00022         worker.finalize()
00023 
00024     except Exception, instance:
00025         sdutil.process_exception(instance)
00026         raise Exception, instance
00027 
00028 class sdplot_worker(sdutil.sdtask_template):
00029     def __init__(self, **kwargs):
00030         super(sdplot_worker,self).__init__(**kwargs)
00031 
00032     def initialize_scan(self):
00033         isScantable = is_scantable(self.infile)
00034 
00035         #load the data without time/pol averaging
00036         sorg = sd.scantable(self.infile,average=self.scanaverage,antenna=self.antenna)
00037 
00038         rfset = (self.restfreq != '') and (self.restfreq != [])
00039         doCopy = (self.frame != '') or (self.doppler != '') or rfset \
00040                  or (self.fluxunit != '' and self.fluxunit != sorg.get_fluxunit()) \
00041                  or (self.specunit != '' and self.specunit != sorg.get_unit())
00042         doCopy = doCopy and isScantable
00043 
00044         # A scantable selection
00045         if type(self.scanlist) == str:
00046             self.scanlist = sorg.parse_idx_selection("SCAN",self.scanlist)
00047         if type(self.iflist) == str:
00048             self.iflist = sorg.parse_idx_selection("IF",self.iflist)
00049         if type(self.pollist) == str:
00050             self.pollist = sorg.parse_idx_selection("POL",self.pollist)
00051         if type(self.beamlist) == str:
00052             self.beamlist = sorg.parse_idx_selection("BEAM",self.beamlist)
00053         sel = self.get_selector()
00054         sorg.set_selection(sel)
00055         self.ssel=sel.__str__()
00056         del sel
00057 
00058         # Copy scantable when usign disk storage not to modify
00059         # the original table.
00060         if doCopy and self.is_disk_storage:
00061             self.scan = sorg.copy()
00062         else:
00063             self.scan = sorg
00064         del sorg
00065 
00066     def execute(self):
00067         # set various attributes to self.scan
00068         self.set_to_scan()
00069 
00070         # Average data if necessary
00071         self.scan = sdutil.doaverage(self.scan, self.scanaverage, self.timeaverage, self.tweight, self.polaverage, self.pweight)
00072 
00073         # Reload plotter if necessary
00074         sd.plotter._assert_plotter(action="reload")
00075 
00076         # Set subplot layout
00077         if self.subplot > 10:
00078             row = int(self.subplot/10)
00079             col = (self.subplot % 10)
00080             sd.plotter.set_layout(rows=row,cols=col,refresh=False)
00081         else:
00082             if self.subplot > -1:
00083                 casalog.post(("Invalid subplot value, %d, is ignored. It should be in between 11 and 99." % self.subplot),priority="WARN")
00084             sd.plotter.set_layout(refresh=False)
00085 
00086         # Set subplot margins
00087         if self.margin != sd.plotter._margins:
00088             sd.plotter.set_margin(margin=self.margin,refresh=False)
00089 
00090         # Actual plotting
00091         getattr(self,'plot_%s'%(self.plottype))()
00092 
00093     def save(self):
00094         # Hardcopy
00095         if (self.outfile != '' ) and not ( self.plottype in ['azel','pointing']):
00096             # currently no way w/o screen display first
00097             sd.plotter.save(self.outfile)
00098         
00099     def plot_pointing(self):
00100         kw = {'scan': self.scan}
00101         if self.outfile != '': kw['outfile'] = self.outfile
00102         sd.plotter.plotpointing(**kw)
00103         
00104         self.__print_header(asaplot=False)
00105 
00106     def plot_azel(self):
00107         kw = {'scan': self.scan}
00108         if self.outfile != '': kw['outfile'] = self.outfile
00109         sd.plotter.plotazel(**kw)
00110 
00111         self.__print_header(asaplot=False)
00112 
00113     def plot_totalpower(self):
00114         sd.plotter.plottp(self.scan)
00115 
00116         self.__print_header(asaplot=False)
00117 
00118     def plot_grid(self):
00119         refresh=False
00120         (nx,ny,cellx,celly,mapcenter) = self.__get_grid_parameters()
00121         self.__dogrid(nx, ny, cellx, celly, mapcenter)
00122 
00123         # Now set fluxunit, specunit, frame, and doppler
00124         sdutil.set_fluxunit(self.scan, self.fluxunit, telescopeparm='fix')
00125         sdutil.set_spectral_unit(self.scan, self.specunit)
00126         sdutil.set_freqframe(self.scan, self.frame)
00127         sdutil.set_doppler(self.scan, self.doppler)
00128 
00129         if self.scan.nchan()==1:
00130             errmsg="Trying to plot the continuum/total power data in 'spectra' mode, please use other plottype options" 
00131             raise Exception,errmsg
00132 
00133         # Smooth the spectrum (if desired)
00134         self.__dosmooth()
00135 
00136         sd.plotter.set_data(self.scan,refresh=refresh)
00137 
00138         # Set colormap, linestyles, and linewidth of plots
00139         self.__setup_plotter()
00140         # Need to specify center and spacing in the epoch and
00141         # unit of DIRECTION column (currently assuming J2000 and rad)
00142         # TODO: need checking epoch and unit of DIRECTION in scantable
00143         crad = None
00144         spacing = None
00145         if self.center != '':
00146             epoch, ra, dec = center.split(' ')
00147             cme = me.direction(epoch, ra, dec)
00148             crad = [qa.convert(cme['m0'],'rad')['value'], \
00149                     qa.convert(cme['m1'],'rad')['value']]
00150         if cellx != '' and celly != '':
00151             spacing = [qa.convert(cellx, 'rad')['value'], \
00152                        qa.convert(celly, 'rad')['value']]
00153         sd.plotter.plotgrid(center=crad,spacing=spacing,rows=ny,cols=nx)
00154 
00155         self.__print_header(asaplot=True)
00156 
00157     def plot_spectra(self):
00158         asaplot=True
00159         if self.scan.nchan()==1:
00160             errmsg="Trying to plot the continuum/total power data in 'spectra' mode, please use other plottype options" 
00161             raise Exception,errmsg
00162 
00163         # Smooth the spectrum (if desired)
00164         self.__dosmooth()
00165 
00166         # Plot final spectrum
00167         # each IF is separate panel, pols stacked
00168         refresh=False
00169 
00170         # Set colormap, linestyles, and linewidth of plots
00171         self.__setup_plotter()
00172         
00173         # The actual plotting
00174         sd.plotter.set_data(self.scan,refresh=refresh)
00175         sd.plotter.plot()
00176 
00177         # Line catalog
00178         dolinc=False
00179         if  self.linecat != 'none' and self.linecat != '':
00180             self.__overlay_linecatalog()
00181 
00182         self.__print_header(asaplot=True)
00183 
00184     def __setup_plotter(self):
00185         refresh = False
00186         sd.plotter.set_mode(stacking=self.stack,panelling=self.panel,refresh=refresh)
00187         ncolor = 0
00188         if self.colormap != 'none': 
00189             colmap = colormap
00190             ncolor=len(colmap.split())
00191         elif self.linestyles == 'none': 
00192             colmap = "green red black cyan magenta orange blue purple yellow pink"
00193             ucm = sd.rcParams['plotter.colours']
00194             if isinstance(ucm,str) and len(ucm) > 0: colmap = ucm
00195             ncolor=len(colmap.split())
00196             del ucm
00197         else: colmap=None
00198 
00199         if self.linestyles != 'none': lines=linestyles
00200         elif ncolor <= 1: 
00201             lines = "line dashed dotted dashdot"
00202             uls = sd.rcParams['plotter.linestyles']
00203             if isinstance(uls,str) and len(uls) > 0: lines = uls
00204             del uls
00205         else: lines=None
00206 
00207         if isinstance(self.linewidth,int) or isinstance (self.linewidth,float):
00208             lwidth = self.linewidth
00209         else:
00210             casalog.post( "Invalid linewidth. linewidth is ignored and set to 1.", priority = 'WARN' )
00211             lwidth = 1
00212 
00213         # set plot colors
00214         if colmap is not None:
00215             if ncolor > 1 and lines is not None:
00216                 casalog.post( "'linestyles' is valid only for single colour plot.\n...Ignoring 'linestyles'.", priority = 'WARN' )
00217             sd.plotter.set_colors(colmap,refresh=refresh)
00218         else:
00219             if lines is not None:
00220                 tmpcol="black"
00221                 #print "INFO: plot colour is set to '",tmpcol,"'"
00222                 casalog.post( "plot colour is set to '"+tmpcol+"'" )
00223                 sd.plotter.set_colors(tmpcol,refresh=refresh)
00224         # set linestyles and/or linewidth
00225         # so far, linestyles can be specified only if a color is assigned
00226         #if lines is not None or linewidth is not None:
00227         #        sd.plotter.set_linestyles(lines, linewidth,refresh=refresh)
00228         sd.plotter.set_linestyles(lines, lwidth,refresh=refresh)
00229         # Plot red x-axis at y=0 (currently disabled)
00230         # sd.plotter.axhline(color='r',linewidth=2)
00231         sd.plotter.set_histogram(hist=self.histogram,refresh=refresh)
00232 
00233         # Set axis ranges (if requested)
00234         if len(self.flrange)==1:
00235             #print "flrange needs 2 limits - ignoring"
00236             casalog.post( "flrange needs 2 limits - ignoring" )
00237         if len(self.sprange)==1:
00238             #print "sprange needs 2 limits - ignoring"
00239             casalog.post( "sprange needs 2 limits - ignoring" )
00240         if ( len(self.sprange) > 1 ):
00241             if ( len(self.flrange) > 1 ):
00242                 sd.plotter.set_range(self.sprange[0],self.sprange[1],self.flrange[0],self.flrange[1],refresh=refresh)
00243             else:
00244                 sd.plotter.set_range(self.sprange[0],self.sprange[1],refresh=refresh)
00245         elif ( len(self.flrange) > 1 ):
00246             sd.plotter.set_range(ystart=self.flrange[0],yend=self.flrange[1],refresh=refresh)
00247         else:
00248             # Set default range explicitly (in case range was ever set)
00249             sd.plotter.set_range(refresh=refresh)
00250 
00251         # legend position
00252         loc=1
00253         if self.plotstyle: loc=legendloc
00254         sd.plotter.set_legend(mode=loc,refresh=refresh) 
00255 
00256     def __overlay_linecatalog(self):
00257         # Use jpl catalog for now (later allow custom catalogs)
00258 
00259         casapath=os.environ['CASAPATH'].split()
00260         catpath=casapath[0]+'/data/catalogs/lines'
00261         catname=catpath+'/jpl_asap.tbl'
00262         # TEMPORARY: hard-wire to my version
00263         # catname='/home/sandrock/smyers/NAUG/Tasks/jpl.tbl'
00264         # FOR LOCAL CATALOGS:
00265         # catname='jpl.tbl'
00266         try:
00267             linc=sd.linecatalog(catname)
00268             dolinc=True
00269         except:
00270             casalog.post( "Could not find catalog at "+catname, priority = "WARN" )
00271             dolinc=False
00272         if ( dolinc ):
00273             if ( len(self.sprange)>1 ):
00274                 if ( specunit=='GHz' or specunit=='MHz' ):
00275                     linc.set_frequency_limits(sprange[0],sprange[1],specunit)
00276                 else:
00277                     casalog.post( "sd.linecatalog.set_frequency_limits accepts onlyGHz and MHz", priority = 'WARN' )
00278                     casalog.post( "continuing without sprange selection on catalog", priority = 'WARN' )
00279             if ( self.linecat != 'all' and self.linecat != 'ALL' ):
00280                 # do some molecule selection
00281                 linc.set_name(self.linecat)
00282             # Plot up the selected part of the line catalog
00283             # use doppler offset
00284             sd.plotter.plot_lines(linc,doppler=self.linedop)
00285             del linc
00286 
00287     def __print_header(self, asaplot=False):
00288         # List observation header
00289         if self.header and (not self.plotstyle or self.margin==[]):
00290             # set margin for the header
00291             sd.plotter._plotter.figure.subplots_adjust(top=0.8)
00292         datname='Data File:     '+self.infile
00293         if self.plottype == 'grid':
00294              datname += " (gridded)"
00295         sd.plotter.print_header(plot=(self.header and asaplot),fontsize=self.headsize,
00296                                 logger=True,selstr=self.ssel,extrastr=datname)        
00297 
00298     def __dosmooth(self):
00299         if self.kernel == '': self.kernel = 'none'
00300         if ( self.kernel != 'none' and (not (self.kwidth<=0 and self.kernel!='hapnning'))):
00301             casalog.post( "Smoothing spectrum with kernel "+self.kernel )
00302             self.scan.smooth(self.kernel,self.kwidth)
00303 
00304     def __dogrid(self, nx, ny, cellx, celly, center):
00305         # Do only the first IFNO and POLNO
00306         ifnos = self.scan.getifnos()
00307         polnos = self.scan.getpolnos()
00308         if len(ifnos) > 1 or len(polnos) > 1:
00309             casalog.post("Only the first IFNO (%d) and POLNO (%d) is plotted." % (ifnos[0], polnos[0]),priority="WARN")
00310 
00311 
00312         gridder = sd.asapgrid2(self.scan)
00313         del self.scan
00314         gridder.setIF(ifnos[0])
00315         gridder.setPolList([polnos[0]])
00316         gridder.defineImage(center=center,cellx=cellx,celly=celly,nx=nx,ny=ny)
00317         gridder.setFunc(func='BOX')
00318         gridder.setWeight('uniform')
00319         gridder.grid()
00320 
00321         self.scan = gridder.getResult()
00322         del gridder
00323 
00324     def __get_grid_parameters(self):
00325         (nx,ny) = self.__get_mapsize()
00326         (cellx,celly) = sdutil.get_cellx_celly(self.cell)
00327         mapcenter = sdutil.get_map_center(self.center)
00328         return (nx,ny,cellx,celly,mapcenter)
00329         
00330     def __get_mapsize(self):
00331         if self.subplot < 11:
00332             casalog.post("Setting default subplot layout (1x1).",priority="WARN")
00333             self.subplot = 11
00334         nx = (self.subplot % 10)
00335         ny = int(self.subplot/10)
00336         return (nx,ny)
00337