diff --git a/Owlcat/Console.py b/Owlcat/Console.py new file mode 100644 index 0000000..4d2e32d --- /dev/null +++ b/Owlcat/Console.py @@ -0,0 +1,43 @@ +# -*- coding: utf-8 -*- + +import sys +import time +import curses + +# Get width of terminal +# I couldn't be bothered to get newline mode working properly when curses is active, +# so I uses curses.wrapper() do init curses, get the width, then close curses. +def get_width (scr): + global _width; + _width = scr.getmaxyx()[1] - 1; +curses.wrapper(get_width); + +def timestamp (time_start,format="%H:%M:%S:"): + return time.strftime(format,time.gmtime(time.time()-time_start)); + +class Reporter (object): + """A Reporter is used to make progress reports on the console, with optional timestamps.""" + def __init__ (self,timestamp=False): + self.time_start = timestamp and time.time(); + + def overprint (self,message): + return self.pprint(message+"\r"); + + def pprint (self,message,newline=True): + # print message + if self.time_start: + message = "%s %s"%(timestamp(self.time_start),message); + if message[-1] == "\r": + endline = "\r"; + else: + endline = "\n\r"; + # cut message at width of terminal, and pad with spaces + sys.stdout.write("%-*.*s"%(_width,_width,message.strip())); + sys.stdout.write(endline); + sys.stdout.flush(); + + def __call__ (self,*args): + self.pprint(" ".join(args)); + + + \ No newline at end of file diff --git a/Owlcat/Flagger.py b/Owlcat/Flagger.py new file mode 100644 index 0000000..50bf578 --- /dev/null +++ b/Owlcat/Flagger.py @@ -0,0 +1,1098 @@ +# -*- coding: utf-8 -*- +import numpy +import Timba.dmi +import re +import tempfile +import os + +import Meow +import Meow.MSUtils +import Purr.Pipe + +from Meow.MSUtils import TABLE + +_gli = Meow.MSUtils.find_exec('glish'); +if _gli: + _GLISH = 'glish'; + Meow.dprint("Calico flagger: found %s, autoflag should be available"%_gli); +else: + _GLISH = None; + Meow.dprint("Calico flagger: glish not found, autoflag will not be available"); + +_addbitflagcol = Meow.MSUtils.find_exec('addbitflagcol'); + +# Various argument-formatting methods to use with the Flagger.AutoFlagger class +# below. These really should be static methods of the class, but that doesn't work +# with Python (specifically, I cannot include them into static member dicts) +def _format_nbins (bins,argname): + if isinstance(bins,(list,tuple)) and len(bins) == 2: + return str(list(bins)); + else: + return "%d"%bins; + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,bins); + +def _format_bool (arg,argname): + if arg: return "T"; + else: return "F"; + +def _format_index (arg,argname): + if isinstance(arg,int): + if arg<0: + return str(arg); + else: + return str(arg+1); + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,arg); + +def _format_list (arg,argname): + if isinstance(arg,str): + return "'%s'"%arg; + elif isinstance(arg,(list,tuple)): + return "[%s]"%(','.join([_format_list(x,argname) for x in arg])); + elif isinstance(arg,bool): + return _format_bool(arg,argname); + elif isinstance(arg,(int,float)): + return str(arg); + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,arg); + +def _format_ilist (arg,argname): + if isinstance(arg,str): + return "'%s'"%arg; + elif isinstance(arg,(list,tuple)): + return "[%s]"%(','.join([_format_ilist(x,argname) for x in arg])); + elif isinstance(arg,bool): + return _format_bool(arg,argname); + elif isinstance(arg,int): + return _format_index(arg,argname); + elif isinstance(arg,float): + return str(arg); + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,arg); + +def _format_plotchan (arg,argname): + if isinstance(arg,bool): + return _format_bool(arg,argname); + else: + return _format_index(arg,argname); + +def _format_2N (arg,argname): + if isinstance(arg,(list,tuple)): + return "%s"%arg; + elif Timba.dmi.is_array(arg): + if arg.dtype == Timba.array.int32: + a = arg + 1; + else: + a = arg; + if a.ndim == 1: + return str(list(a)); + elif a.ndim == 2 and a.shape[1] == 2: + return "[%s]"%(','.join([ str(list(a[i,:])) for i in range(a.shape[0]) ])); + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,arg); + +def _format_clip (arg,argname): + if isinstance(arg,(list,tuple)): + if not isinstance(arg[0],(list,tuple)): + arg = [ arg ]; + recfields = []; + for i,triplet in enumerate(arg): + if not isinstance(triplet,(list,tuple)) or len(triplet) != 3: + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,arg); + expr,minval,maxval = triplet; + subfields = [ "expr='%s'"%expr ]; + if minval is not None: + subfields.append("min=%g"%minval); + if maxval is not None: + subfields.append("max=%g"%maxval); + recfields.append("a%d=[%s]"%(i,','.join(subfields))); + return "[%s]"%','.join(recfields); + raise TypeError,"invalid value for '%s' keyword (%s)"%(argname,arg); + +class Flagger (Timba.dmi.verbosity): + def __init__ (self,msname,verbose=0,chunksize=200000): + Timba.dmi.verbosity.__init__(self,name="Flagger"); + self.set_verbose(verbose); + if not TABLE: + raise RuntimeError,"No tables module found. Please install pyrap and casacore"; + self.msname = msname; + self.ms = None; + self.readwrite = False; + self.chunksize = chunksize; + self._reopen(); + + def close (self): + if self.ms: + self.dprint(2,"closing MS",self.msname); + self.ms.close(); + self.ms = self.flagsets = None; + + def _reopen (self,readwrite=False): + if self.ms is None: + self.ms = ms = TABLE(self.msname,readonly=not readwrite); + self.readwrite = readwrite; + self.dprintf(1,"opened MS %s, %d rows\n",ms.name(),ms.nrows()); + self.has_bitflags = 'BITFLAG' in ms.colnames(); + self.flagsets = Meow.MSUtils.get_flagsets(ms); + self.dprintf(1,"flagsets are %s\n",self.flagsets.names()); + self.purrpipe = Purr.Pipe.Pipe(self.msname); + elif self.readwrite != readwrite: + self.ms = TABLE(self.msname,readonly=not readwrite); + self.readwrite = readwrite; + return self.ms; + + def add_bitflags (self,wait=True,purr=True): + if not self.has_bitflags: + global _addbitflagcol; + if not _addbitflagcol: + raise RuntimeError,"cannot find addbitflagcol utility in $PATH"; + self.close(); + self.dprintf(1,"running addbitflagcol\n"); + if os.spawnvp(os.P_WAIT,_addbitflagcol,['addbitflagcol',self.msname]): + raise RuntimeError,"addbitflagcol failed"; + # report to purr + if purr: + self.purrpipe.title("Flagging").comment("Adding bitflag columns."); + # reopen and close to pick up new BITFLAG column + self._reopen(); + self.close(); + + def remove_flagset (self,*fsnames,**kw): + self._reopen(True); + mask = self.flagsets.remove_flagset(*fsnames); + # report to purr + if kw.get('purr',True): + self.purrpipe.title("Flagging").comment("Removing flagset(s) %s."%(','.join(fsnames))); + self.unflag(mask,purr=False); + + def flag (self,flag=1,**kw): + kw = kw.copy(); + if kw.setdefault('purr',True): + if isinstance(flag,int): + fset = "bitflag %x"%flag; + else: + fset = "flagset %s"%flag; + self.purrpipe.title("Flagging").comment("Flagging %s"%fset,endline=False); + kw['flag'] = flag; + kw['unflag'] = 0; + return self._flag(**kw); + + def unflag (self,unflag=-1,**kw): + kw = kw.copy(); + if kw.setdefault('purr',True): + if isinstance(unflag,int): + fset = "bitflag %x"%unflag; + else: + fset = "flagset %s"%unflag; + self.purrpipe.title("Flagging").comment("Unflagging %s"%fset,endline=False); + kw['flag'] = 0; + kw['unflag'] = unflag; + return self._flag(**kw); + + def transfer (self,flag=1,replace=False,*args,**kw): + unflag = (replace and flag) or 0; + kw.setdefault('purr',True); + if kw['purr']: + if isinstance(flag,int): + fset = "bitflag %x"%flag; + else: + fset = "flagset %s"%flag; + self.purrpipe.title("Flagging").comment("Transferring FLAG/FLAG_ROW to %s"%fset,endline=False); + return self._flag(flag=flag,unflag=unflag,transfer=True,*args,**kw); + + def get_stats(self,flag=0,legacy=False,**kw): + kw.setdefault('purr',True); + if kw['purr']: + if isinstance(flag,int) and flag: + fset = "bitflag %x"%flag; + else: + fset = "flagset %s"%flag; + if legacy: + fset += ", plus FLAG/FLAG_ROW"; + self.purrpipe.title("Flagging").comment("Getting flag stats for %s"%fset,endline=False); + stats = self._flag(flag=flag,get_stats=True,include_legacy_stats=legacy,**kw); + msg = "%.2f%% of rows and %.2f%% of correlations are flagged."%(stats[0]*100,stats[1]*100); + if kw['purr']: + self.purrpipe.comment(msg); + self.dprint(1,"stats: ",msg); + return stats; + + def _get_bitflag_col (self,ms,row0,nrows,shape=None): + """helper method. Gets the bitflag column at the specified location. On error (presumably, + column is missing), returns zero array of specified shape. If shape is not specified, + queries the DATA column to obtain it."""; + try: + return ms.getcol('BITFLAG',row0,nrows); + except: + if not shape: + shape = ms.getcol('DATA',row0,nrows).shape; + return numpy.zeros(shape,dtype=numpy.int32); + + def _get_submss (self,ms,ddids=None): + """Helper method. Splits MS into subsets by DATA_DESC_ID. + Returns list of (ddid,nrows,subms) tuples, where subms is a subset of the MS with the given DDID, + and nrows is the number of rows in preceding submss (which is needed for progress stats) + """; + # build up list of (ddid,nrows,subms) pairs. + # subms is a subset of the MS with the given DDID + # nrows is the number of rows in preceding submss (need for progress stats) + sub_mss = []; + nrows = 0; + if ddids is None: + ddids = range(TABLE(ms.getkeyword('DATA_DESCRIPTION')).nrows()); + for ddid in ddids: + subms = ms.query("DATA_DESC_ID==%d"%ddid); + sub_mss.append((ddid,nrows,subms)); + nrows += subms.nrows(); + return sub_mss; + + def _flag (self, + flag=1, # set this flagmask (or flagset name) or + unflag=0, # clear this flagmask (or flagset name) + create=False, # if True and 'flag' is a string, creates new flagset as needed + fill_legacy=None, # if not None, legacy flags will be filled using the specified flagmask + transfer=False, # if True: sets transfer mode. In transfer mode, legacy flags (within + # the specified subset) are transferred to the given bitflag + get_stats=False, # if True: sets "get stats" mode. Instead of flagging, counts and + # returns a count of raised flags (matching the flagmask) in the subset + include_legacy_stats=False, # if True and get_stats=True, includes legacy flags in the stats + # the following options restrict the subset to be flagged. Effect is cumulative. + ddid=None,fieldid=None, # DATA_DESC_ID or FIELD_ID, single int or list + channels=None, # channel subset (index, list of indices, or slice object) + multichan=None, # list of channel subsets (as an alternative to specifying just one) + corrs=None, # correlation subset (index, list of indices, or slice object) + antennas=None, # list of antennas + baselines=None, # list of baselines (as ip,iq pairs) + time=None, # time range as (min,max) pair (can use None for either min or max) + reltime=None, # relative time (rel. to start of MS) range as (min,max) pair + taql=None, # additional TaQL string + clip_above=None, # restrict flagged subset to abs(data)>clip_above + clip_below=None, # and abs(data)=%g"%t0); + if t1 is not None: + queries.append("TIME<=%g"%t1); + if reltime is not None: + t0,t1 = reltime; + time0 = self.ms.getcol('TIME',0,1)[0]; + if t0 is not None: + queries.append("TIME>=%f"%(time0+t0)); + if t1 is not None: + queries.append("TIME<=%f"%(time0+t1)); + + # form up TaQL string, and extract subset of table + if queries: + query = "( " + " ) && ( ".join(queries)+" )"; + purr and self.purrpipe.comment("; effective MS selection is \"%s\""%query,endline=False); + self.dprintf(2,"selection string is %s\n",query); + ms = ms.query(query); + self.dprintf(2,"query reduces MS to %d rows\n",ms.nrows()); + else: + self.dprintf(2,"no selection applied\n"); + # check list of baselines + if baselines: + baselines = [ (int(p),int(q)) for p,q in baselines ]; + purr and self.purrpipe.comment("; baseline subset is %s"% + " ".join(["%d-%d"%(p,q) for p,q in baselines]), + endline=False); + + # This will be true if only whole rows are being selected for. If channel/correlation/clipping + # criteria are supplied, this will be set to False below + clip = not ( clip_above is None and clip_below is None and + clip_fm_above is None and clip_fm_below is None); + flagrows = not clip; + # form up channel and correlation slicers + # multichan may specify multiple channel subsets. If not specified, + # then channels specifies a single subset. In any case, in the end multichan + # will contain a list of the current channel selection + if multichan is None: + if channels is None: + multichan = [ numpy.s_[:] ]; + else: + flagrows = False; + multichan = [ channels ]; + purr and self.purrpipe.comment("; channels are %s"%channels,endline=False); + else: + purr and self.purrpipe.comment("; channels are %s"%(', '.join(map(str,multichan))),endline=False); + flagrows = False; + self.dprintf(2,"channel selection is %s\n",multichan); + if corrs is None: + corrs = numpy.s_[:]; + else: + purr and self.purrpipe.comment("; correlations %s"%corrs,endline=False); + flagrows = False; + # putt comment into purrpipe + purr and self.purrpipe.comment("."); + self.dprintf(2,"correlation selection is %s\n",corrs); + stat_rows_nfl = stat_rows = stat_pixels = stat_pixels_nfl = 0; + # make list of sub-MSs by DDID + sub_mss = self._get_submss(ms,ddids); + nrow_tot = ms.nrows(); + # go through rows of the MS in chunks + for ddid,irow_prev,ms in sub_mss: + self.dprintf(2,"processing MS subset for ddid %d\n",ddid); + if progress_callback: + progress_callback(irow_prev,nrow_tot); + for row0 in range(0,ms.nrows(),self.chunksize): + if progress_callback: + progress_callback(irow_prev+row0,nrow_tot); + nrows = min(self.chunksize,ms.nrows()-row0); + self.dprintf(2,"flagging rows %d:%d\n",row0,row0+nrows-1); + # get mask of matching baselines + if baselines: + # init mask of all-false + rowmask = numpy.zeros(nrows,dtype=numpy.bool); + a1 = ms.getcol('ANTENNA1',row0,nrows); + a2 = ms.getcol('ANTENNA2',row0,nrows); + # update mask + for p,q in baselines: + rowmask |= (a1==p) & (a2==q); + self.dprintf(2,"baseline selection leaves %d rows\n",len(rowmask.nonzero()[0])); + # else select all rows + else: + rowmask = numpy.s_[:]; + self.dprintf(2,"no baseline selection applied, flagging %d rows\n",nrows); + # form up subsets for channel/correlation selector + subsets = [ (rowmask,ch,corrs) for ch in multichan ]; + # first, handle statistics mode + if get_stats: + # collect row stats + if include_legacy_stats: + lfr = ms.getcol('FLAG_ROW',row0,nrows)[rowmask]; + lf = ms.getcol('FLAG',row0,nrows); + else: + lfr = lf = 0; + if flag: + bfr = ms.getcol('BITFLAG_ROW',row0,nrows)[rowmask]; + lfr = lfr + ((bfr&flag)!=0); + bf = self._get_bitflag_col(ms,row0,nrows); + # size seems to be a method or an attribute depending on numpy version :( + stat_rows += (callable(lfr.size) and lfr.size()) or lfr.size; + stat_rows_nfl += lfr.sum(); + for subset in subsets: + if include_legacy_stats: + lfm = lf[subset]; + else: + lfm = 0; + if flag: + lfm = lfm + (bf[subset]&flag)!=0; + # size seems to be a method or an attribute depending on numpy version :( + stat_pixels += (callable(lfm.size) and lfm.size()) or lfm.size; + stat_pixels_nfl += lfm.sum(); + # second, handle transfer-flags mode + elif transfer: + bf = ms.getcol('BITFLAG_ROW',row0,nrows); + bfm = bf[rowmask]; + if unflag: + bfm &= ~unflag; + lf = ms.getcol('FLAG_ROW',row0,nrows)[rowmask]; + bf[rowmask] = numpy.where(lf,bfm|flag,bfm); + # size seems to be a method or an attribute depending on numpy version :( + stat_rows += (callable(lf.size) and lf.size()) or lf.size; + stat_rows_nfl += lf.sum(); + ms.putcol('BITFLAG_ROW',bf,row0,nrows); + lf = ms.getcol('FLAG',row0,nrows); + bf = self._get_bitflag_col(ms,row0,nrows,lf.shape); + for subset in subsets: + bfm = bf[subset]; + if unflag: + bfm &= ~unflag; + lfm = lf[subset] + bf[subset] = numpy.where(lfm,bfm|flag,bfm); + # size seems to be a method or an attribute depending on numpy version :( + stat_pixels += (callable(lfm.size) and lfm.size()) or lfm.size; + stat_pixels_nfl += lfm.sum(); + ms.putcol('BITFLAG',bf,row0,nrows); + # else, are we flagging whole rows? + elif flagrows: + if self.has_bitflags: + bfr = ms.getcol('BITFLAG_ROW',row0,nrows); + bf = self._get_bitflag_col(ms,row0,nrows); + if unflag: + bfr[rowmask] &= ~unflag; + bf[rowmask,:,:] &= ~unflag; + if flag: + bfr[rowmask] |= flag; + bf[rowmask,:,:] |= flag; + ms.putcol('BITFLAG_ROW',bfr,row0,nrows); + ms.putcol('BITFLAG',bf,row0,nrows); + if fill_legacy is not None: + lfr = ms.getcol('FLAG_ROW',row0,nrows); + lf = ms.getcol('FLAG',row0,nrows); + lfr[rowmask] = ( (bfr[rowmask]&fill_legacy) !=0 ); + lf[rowmask,:,:] = ( (bf[rowmask]&fill_legacy) !=0 ); + ms.putcol('FLAG_ROW',lfr,row0,nrows); + ms.putcol('FLAG',lf,row0,nrows); + else: + lfr = ms.getcol('FLAG_ROW',row0,nrows); + lf = ms.getcol('FLAG',row0,nrows); + lfr[rowmask] = (flag!=0); + lf[rowmask,:,:] = (flag!=0); + ms.putcol('FLAG_ROW',lfr,row0,nrows); + ms.putcol('FLAG',lf,row0,nrows); + # else flagging individual correlations or channels + else: + # get flags (for clipping purposes) + lf = ms.getcol('FLAG',row0,nrows); + # 'mask' is what needs to be flagged/unflagged. Start with empty mask. + mask = numpy.zeros(lf.shape,bool); + # then fill in subsets + for subset in subsets: + mask[subset] = True; + # get clipping mask, if amplitude clipping is in effect + if clip: + datacol = ms.getcol(clip_column,row0,nrows); + clip_mask = numpy.ones(datacol.shape,bool); + if clip_above is not None: + clip_mask &= abs(datacol)>clip_above; + if clip_below is not None: + clip_mask &= abs(datacol) 1: + # mask data column with subset, and with legacy flags + datacol = numpy.ma.masked_array(abs(datacol),(~mask)|lf,fill_value=0).mean(1); + if clip_fm_above is not None: + clip_mask &= (datacol>clip_fm_above)[:,numpy.newaxis,...]; + if clip_fm_below is not None: + clip_mask &= (datacolX + data_below=None, # and abs(data)=%g"%t0); + if t1 is not None: + queries.append("TIME<=%g"%t1); + if reltime is not None: + t0,t1 = reltime; + time0 = self.ms.getcol('TIME',0,1)[0]; + if t0 is not None: + queries.append("TIME>=%f"%(time0+t0)); + if t1 is not None: + queries.append("TIME<=%f"%(time0+t1)); + # form up TaQL string, and extract subset of table + if queries: + query = "( " + " ) && ( ".join(queries)+" )"; + purr and self.purrpipe.comment("; effective MS selection is \"%s\""%query,endline=False); + self.dprintf(2,"selection string is %s\n",query); + ms = ms.query(query); + self.dprintf(2,"query reduces MS to %d rows\n",ms.nrows()); + else: + self.dprintf(2,"no selection applied\n"); + # check list of baselines + if baselines: + baselines = [ (int(p),int(q)) for p,q in baselines ]; + purr and self.purrpipe.comment("; baseline subset is %s"% + " ".join(["%d-%d"%(p,q) for p,q in baselines]), + endline=False); + # helper func to parse the channels/corrs/timeslots arguments + def make_slice_list (selection,parm): + if not selection: + return [ numpy.s_[:] ]; + if isinstance(selection,(int,slice)): + return make_slice_list([selection],parm); + if not isinstance(selection,(list,tuple)): + raise TypeError,"invalid %s selection: %s"%(parm,selection); + sellist = []; + for sel in selection: + if isinstance(sel,int): + sellist.append(slice(sel,sel+1)); + elif isinstance(sel,slice): + sellist.append(sel); + else: + raise TypeError,"invalid %s selection: %s"%(parm,selection); + return sellist; + # parse the arguments + channels = make_slice_list(channels,'channels'); + corrs = make_slice_list(corrs,'corrs'); + purr and self.purrpipe.comment("; channels are %s"%channels,endline=False); + purr and self.purrpipe.comment("; correlations are %s"%corrs,endline=False); + # put comment into purrpipe + purr and self.purrpipe.comment("."); + # + flagsubsets = flagmask is not None or flagmask_all is not None or flagmask_none is not None; + dataclip = data_above is not None or data_below is not None or \ + data_fm_above is not None or data_fm_below is not None; + # make list of sub-MSs by DDID + sub_mss = self._get_submss(ms,ddids); + nrow_tot = ms.nrows(); + # go through rows of the MS in chunks + for ddid,irow_prev,ms in sub_mss: + self.dprintf(2,"processing MS subset for ddid %d\n",ddid); + if progress_callback: + progress_callback(irow_prev,nrow_tot); + for row0 in range(0,ms.nrows(),self.chunksize): + if progress_callback: + progress_callback(irow_prev+row0,nrow_tot); + nrows = min(self.chunksize,ms.nrows()-row0); + self.dprintf(2,"processing rows %d:%d (%d rows total)\n",row0,row0+nrows-1,nrows); + # apply baseline selection to the mask + if baselines: + # rowmask will be True for all selected rows + rowmask = numpy.zeros(nrows,bool); + a1 = ms.getcol('ANTENNA1',row0,nrows); + a2 = ms.getcol('ANTENNA2',row0,nrows); + # update mask + for p,q in baselines: + rowmask |= (a1==p) & (a2==q); + self.dprintf(2,"baseline selection leaves %d rows\n",nr); + # else select all rows + else: + # rowmask will be True for all selected rows + rowmask = numpy.ones(nrows,bool); + # read legacy flags to get a datashape + lf = ms.getcol('FLAG',row0,nrows); + datashape = lf.shape; + nv_per_row = datashape[1]*datashape[2]; + # rowflags and visflags will be constructed on-demand below. Make helper functions for this + self._rowflags = self._visflags = None; + def rowflags (): + if self._rowflags is None: + # read legacy flags and convert them to bitmask, then add bitflags + lfr = ms.getcol('FLAG_ROW',row0,nrows); + self._rowflags = lfr*self.LEGACY; + if self.has_bitflags: + self._rowflags |= ms.getcol('BITFLAG_ROW',row0,nrows); + return self._rowflags; + def visflags (): + if self._visflags is None: + self._visflags = lf*self.LEGACY; + if self.has_bitflags: + bf = ms.getcol('BITFLAG',row0,nrows); + self._bitflag_dtype = bf.dtype; + self._visflags |= bf; + return self._visflags; + # apply stats + nr = rowmask.sum(); + sel_nrow += nr; + nv = nr*nv_per_row; + sel_nvis += nv; + self.dprintf(2,"Row subset (data selection) leaves %d rows and %d visibilities\n",nr,nv); + # get subset C + # vismask will be True for all selected visibilities + vismask = numpy.zeros(datashape,bool); + for channel_slice in channels: + for corr_slice in corrs: + vismask[rowmask,channel_slice,corr_slice] = True; + nv = vismask.sum(); + nvis_A += vismask.sum(); + self.dprintf(2,"subset A (freq/corr slicing) leaves %d visibilities\n",nv); + # read flags if selecting subset D on them (and also if clipping data) + if flagsubsets: + vf = visflags(); + # apply them to the rowmask + if flagmask is not None: + vismask &= ( (vf&flagmask) != 0 ); + if flagmask_all is not None: + vismask &= ( (vf&flagmask_all) == flagmask_all ); + if flagmask_none is not None: + vismask &= ( (vf&flagmask_none) == 0 ); + nv = vismask.sum(); + nvis_B += nv; + self.dprintf(2,"subset B (flag-based selection) leaves %d visibilities\n",nv); + # now apply clipping + if dataclip: + datacol = ms.getcol(data_column,row0,nrows); + # make it a masked array: mask out stuff not in vismask + datamask = ~vismask; + # and mask stuff in data_flagmask + if data_flagmask is not None: + datamask |= ( (visflags()&data_flagmask)!=0 ); + datacol = numpy.ma.masked_array(datacol,datamask); + # clip on amplitudes + if data_above is not None: + vismask &= abs(datacol)>data_above; + if data_below is not None: + vismask &= abs(datacol)data_fm_above)[:,numpy.newaxis,...]; + if data_fm_below is not None: + vismask &= (datacol Wrote",save; + return fig; + diff --git a/Owlcat/__init__.py b/Owlcat/__init__.py new file mode 100644 index 0000000..eeb0adb --- /dev/null +++ b/Owlcat/__init__.py @@ -0,0 +1,106 @@ +# -*- coding: utf-8 -*- +# +#% $Id: __init__.py 6998 2009-06-26 08:51:15Z cwilliams $ +# +# +# Copyright (C) 2002-2007 +# The MeqTree Foundation & +# ASTRON (Netherlands Foundation for Research in Astronomy) +# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see , +# or write to the Free Software Foundation, Inc., +# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# + +import __main__ +setattr(__main__,"_meow_verbosity",0); + +# import table class +try: + from pyrap_tables import table,tablecopy,tableexists,tabledelete +except: + try: + from pyrap.tables import table,tablecopy,tableexists,tabledelete + except: + print "Failed to import pyrap_tables or pyrap.tables. Please install the pyrap " + "package from http://code.google.com/p/pyrap/, or from a MeqTrees binary repository " + "(see http://www.astron.nl/meqwiki/Downloading)" + raise; + +# pyrap likes to display a lot of these, so shut them off +import warnings +warnings.simplefilter('ignore',DeprecationWarning); + +# list of optional packages which will be added to the include path +_Packages = [ "Cattery" ]; + +# list of locations where packages will be searched for +_PackageLocations = [ "~","~/Frameworks/", + "/usr/local/MeqTrees","/usr/local/lib/MeqTrees","/usr/lib/MeqTrees","/usr/lib64/MeqTrees","/usr/lib32/MeqTrees" + "/usr/local/meqtrees","/usr/local/lib/meqtrees","/usr/lib/meqtrees","/usr/lib64/meqtrees","/usr/lib32/meqtrees" + ]; + +# mapping of package: path. Filled in as we find packages +_packages = {}; + +import sys +import os +import os.path + +def packages (): + """Returns mapping of available packages to their paths"""; + return _packages; + # print "Using %s, set the %s_PATH environment variable to override this."%(path,package.upper()); + +def _tryPackageDir (path,package): + """Tests if path refers to a valid directory, adds it to system include path if so. + Marks package as having this path."""; + if os.path.isdir(path): + sys.path.insert(0,path); + # check for version info + try: + version = ' '.join(file(os.path.join(path,'version_info'))); + except: + version = 'no version info'; + # insert into packages + global _packages; + _packages[package] = path,version; + return True; + return False; + +def _setPackagePath (package): + """Finds the given package, by first looking in $MEQTREES_PACKAGE_PATH, then checking for + subdirectories of the standard _PackageLocations list."""; + # check for explicit MEQTREES_PACKAGE_PATH first + varname = 'MEQTREES_%s_PATH'%package.upper() + path = os.environ.get(varname,None); + if path: + if not _tryPackageDir(path,package): + print "Warning: your %s environment variable is set to"%varname; + print "%s, but this is not a valid directory."%path; + print "The %s package will not be available."%package; + return; + # else look in standard places + for path in _PackageLocations: + path = os.path.expanduser(path); + if _tryPackageDir(os.path.join(path,package),package): + return; + # none found + print "Warning: No %s package found."%package; + print "If you have %s in a non-standard location, please set the %s environment"%(package,varname); + print "variable to point to it." + +for pkg in _Packages: + _setPackagePath(pkg); diff --git a/downweigh-redundant-baselines.py b/downweigh-redundant-baselines.py new file mode 100755 index 0000000..57ea9e6 --- /dev/null +++ b/downweigh-redundant-baselines.py @@ -0,0 +1,124 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import sys +import os.path +import re +import time +import math +import Owlcat + +ROWCHUNK = 50000; + +time_start = time.time(); + +try: + from pyrap_tables import table,tablecopy,tableexists,tabledelete +except: + from pyrap.tables import table,tablecopy,tableexists,tabledelete + +def timestamp (format="%H:%M:%S:"): + return time.strftime(format,time.gmtime(time.time()-time_start)); + +def progress (message,newline=True): + sys.stdout.write("%s %-60s%s"%(timestamp(),message,'\n' if newline else '\r')); + sys.stdout.flush(); + +if __name__ == "__main__": + + # setup some standard command-line option parsing + # + from optparse import OptionParser + parser = OptionParser(usage="""%prog: [options] MS""", + description="""Finds redundant baselines and modifies the IMAGING_WEIGHT column of the MS, +reweighting the redundant baselines by 1/n. This compensates for the lack of this option +in the casa imager. For proper treatment of redundant baselines in e.g. WSRT, first make +an image using radial weighting (so that the imager fills in the IMAGING_WEIGHT column +appropriately), then run this script to modify the weights, then run the imager again, but +this time tell it to use the default weights."""); + #parser.add_option("-o","--output",dest="output",type="string", + #help="name of output FITS file"); + #parser.add_option("-z","--zoom",dest="output_quad",type="string", + #help="name of zoomed output FITS file"); + parser.add_option("-t","--tolerance",dest="tolerance",type="float", + help="How close (in meters) two baselines need to be to each other to be considered redundant (default .1)"); + parser.add_option("-I","--ifrs",dest="ifrs",type="string", + help="subset of interferometers to use."); + parser.add_option("-s","--select",dest="select",action="store", + help="additional TaQL selection string. Note that redundant baselines are counted only within the subset " + "given by the --ifrs and --select options."); + parser.set_defaults(tolerance=.1,select="",ifrs=""); + + (options,msnames) = parser.parse_args(); + + if len(msnames) != 1: + parser.error("MS name not supplied."); + + msname = msnames[0]; + ms = table(msname,readonly=False); + + taqls = []; + # get IFR set + import Meow.IfrSet + ifrset = Meow.IfrSet.from_ms(ms); + if options.ifrs: + ifrset = ifrset.subset(options.ifrs); + taqls.append(ifrset.taql_string()); + + if options.select: + taqls.append(options.select); + + if taqls: + select = "( " + " ) &&( ".join(taqls) + " )"; + progress("Applying TaQL selection %s"%select,newline=True); + ms = ms.query(select); + progress("Looking for redundant baselines",newline=True); + ant1 = ms.getcol('ANTENNA1'); + ant2 = ms.getcol('ANTENNA2'); + + IFRS = sorted(set([ (p,q) for p,q in zip(ant1,ant2) ])); + print "%d baselines"%len(IFRS); + groups = []; + + for i,(p,q) in enumerate(IFRS): + bl = ifrset.baseline_vector(p,q); + # see if this baseline is within the tolerance of a previous group's baseline + for ig,(bl0,members) in enumerate(groups): + if abs(bl-bl0).max() < options.tolerance: + members.append((p,q)); + break; + # if none, start a new group + else: + members = [(p,q)]; + ig = len(groups); + groups.append([]); + # update baseline (as mean baseline of group) + length = reduce(lambda x,y:x+y,[ifrset.baseline_vector(*ifr) for ifr in members])/len(members); + groups[ig] = (length,members); + + # make a dictionary of per-IFR weights + have_redundancy = False; + weight = dict([((p,q),1.) for p,q in IFRS]); + for baseline,members in sorted(groups,lambda a,b:cmp((a[0]**2).sum(),(b[0]**2).sum())): + if len(members) > 1: + print "Baseline %dm, %d ifrs: %s"%(round(math.sqrt((baseline**2).sum())),len(members), + " ".join(["%d-%d"%(p,q) for p,q in members])); + have_redundancy = True; + for p,q in members: + weight[p,q] = 1.0/len(members); + if not have_redundancy: + print "No redundant baselines found, nothing to do." + sys.exit(0); + + # apply weights + nrows = ms.nrows(); + for row0 in range(0,nrows,ROWCHUNK): + progress("Applying weights, row %d of %d"%(row0,nrows),newline=False); + row1 = min(nrows,row0+ROWCHUNK); + w = ms.getcol('IMAGING_WEIGHT',row0,row1-row0); + for i in range(row0,row1): + w[i-row0,:] *= weight[ant1[i],ant2[i]]; + ms.putcol('IMAGING_WEIGHT',w,row0,row1-row0); + + progress("Closing MS."); + ms.close(); diff --git a/fitstool.py b/fitstool.py new file mode 100755 index 0000000..fb02dec --- /dev/null +++ b/fitstool.py @@ -0,0 +1,98 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import sys +import pyfits +import re +import os.path + +if __name__ == "__main__": + + # setup some standard command-line option parsing + # + from optparse import OptionParser + parser = OptionParser(usage="""%prog: [options] """); + parser.add_option("-o","--output",dest="output",type="string", + help="name of output FITS file"); + parser.add_option("-f","--force",dest="force",action="store_true", + help="overwrite output file even if it exists"); + parser.add_option("-m","--mean",dest="mean",action="store_true", + help="take mean of input images"); + parser.add_option("-d","--diff",dest="diff",action="store_true", + help="take difference of 2 input images"); + parser.add_option("-z","--zoom",dest="zoom",type="int",metavar="NPIX", + help="zoom into central region of NPIX x NPIX size"); + parser.add_option("-s","--scale",dest="scale",type="float", + help="rescale image values"); + parser.set_defaults(output="",mean=False,zoom=0,scale=1); + + (options,imagenames) = parser.parse_args(); + + if not imagenames: + parser.error("No images specified. Use '-h' for help."); + + print "%d input image(s): %s"%(len(imagenames),", ".join(imagenames)); + images = [ pyfits.open(img) for img in imagenames ]; + updated = False; + + outname = options.output; + autoname = not outname; + if autoname: + outname = re.split('[_]',imagenames[0],1)[-1]; + + if options.diff and options.mean: + parser.error("Cannot do both --mean and --diff at once."); + + if options.diff: + if len(images) != 2: + parser.error("The --diff option requires exactly two input images."); + if autoname: + outname = "diff_" + outname; + print "Computing difference"; + data = images[0][0].data; + data -= images[1][0].data; + updated = True; + + if options.mean: + if autoname: + outname = "mean%d_"%len(images) + outname; + print "Computing mean"; + data = images[0][0].data; + for img in images[1:]: + data += img[0].data; + data /= len(images); + images = [ images[0] ]; + updated = True; + + if options.zoom: + z = options.zoom; + if autoname: + outname = "zoom%d_"%z + outname; + if len(images) > 1: + "Too many input images specified for this operation, at most 1 expected"; + sys.exit(2); + data = images[0][0].data; + nx = data.shape[2]; + ny = data.shape[3]; + zdata = data[:,:,(nx-z)/2:(nx+z)/2,(ny-z)/2:(ny+z)/2]; + print "Making zoomed image of shape","x".join(map(str,zdata.shape)); + images = [ pyfits.PrimaryHDU(zdata) ]; + updated = True; + + if options.scale != 1: + if autoname and not updated: + outname = "rescale_" + outname; + if len(images) > 1: + "Too many input images specified for this operation, at most 1 expected"; + sys.exit(2); + print "Applying scaling factor of %f to image values"%options.scale; + images[0][0].data *= options.scale; + updated = True; + + if updated: + print "Writing output image",outname; + if os.path.exists(outname) and not options.force: + parser.error("Output image exists, rerun with the -f switch to overwrite."); + images[0].writeto(outname,clobber=True); + else: + print "No operations specified. Use --help for help." diff --git a/flag-ms.py b/flag-ms.py new file mode 100755 index 0000000..dd112b6 --- /dev/null +++ b/flag-ms.py @@ -0,0 +1,349 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import os.path +import math +import sys +import re +import traceback + +flagger = parser = ms = msname = None; + +def error (message): + print "%s: %s"%(os.path.basename(sys.argv[0]),message); + sys.exit(1); + +def get_ms (): + global ms; + global msname; + if not ms: + ms = Owlcat.table(msname); + return ms; + +def parse_subset_options (options): + global parser; + global flagger; + global msname; + global ms; + subset = {}; + from Owlcat import Parsing + + # DDID and FIELD_ID + if options.ddid is not None: + try: + subset['ddid'] = map(int,options.ddid.split(",")); + print " ===> DATA_DESC_ID:",subset['ddid']; + except: + parser.error("Invalid -D/--ddid option"); + if options.field is not None: + try: + subset['fieldid'] = map(int,options.field.split(",")); + print " ===> FIELD_ID:",subset['fieldid']; + except: + parser.error("Invalid -F/--field option"); + # taql + taqls = []; + if options.taql: + taqls.append(options.taql); + print " ===> TaQL query:",options.taql; + # channels + if options.channels: + subset['channels'] = map(Parsing.parse_slice,options.channels.split(",")); + print " ===> channels:",subset['channels']; + # corr list + if options.corrs is not None: + try: + subset['corrs'] = map(int,options.corrs.split(',')); + print " ===> correlations:",subset['corrs']; + except: + parser.error("Invalid -X/--corrs option"); + # station list + if options.stations is not None: + try: + subset['antennas'] = map(int,options.stations.split(',')); + print " ===> stations:",subset['antennas']; + except: + parser.error("Invalid -S/--stations option"); + # IFR set + if options.ifrs is not None: + import Meow.IfrSet + ifrset = Meow.IfrSet.from_ms(get_ms()); + # print help and exit + if options.ifrs == "help": + # print help string, but trim away RTF tags + print re.sub("<[^>]+>","",ifrset.subset_doc).replace("<","<").replace(">",">"); + sys.exit(0); + try: + ifrset = ifrset.subset(options.ifrs); + print " ===> ifrs:"," ".join([ifrset.ifr_label(ip,iq) for (ip,p),(iq,q) in ifrset.ifr_index()]); + except: + parser.error("Invalid -I/--ifrs option"); + taqls.append(ifrset.taql_string()); + # clipping + subset['data_column'] = options.data_column; + if options.above is not None: + subset['data_above'] = options.above; + print " ===> select |%s|>%f"%(options.data_column,options.above); + if options.below is not None: + subset['data_below'] = options.below; + print " ===> select |%s|<%f"%(options.data_column,options.below); + if options.fm_above is not None: + subset['data_fm_above'] = options.fm_above; + print " ===> select mean|%s|>%f"%(options.data_column,options.fm_above); + if options.fm_below is not None: + subset['data_fm_below'] = options.fm_below; + print " ===> select mean|%s|<%f"%(options.data_column,options.fm_below); + # join taql queries + if taqls: + subset['taql'] = "( " + " ) && ( ".join(taqls) + " )"; + # fill flag args + for opt in 'data_flagmask','flagmask','flagmask_all','flagmask_none': + subset[opt] = getattr(options,opt); + return subset; + + +if __name__ == "__main__": + + # setup some standard command-line option parsing + # + from optparse import OptionParser,OptionGroup + parser = OptionParser(usage="""%prog: [actions] [options] MS""", + description="Manipulates flags (bitflags and legacy FLAG/FLAG_ROW columns) in the MS. " + "Use the selection options to narrow down a subset of the data, and use the action options " + "to change flags within that subset. Without any action options, statistics on the current " + "selection are printed -- this is useful as a preview of your intended action." + ); + + group = OptionGroup(parser,"Selection by subset"); + group.add_option("-L","--channels",type="string", + help="channel selection, as or ::[:]. 'end' is inclusive. Can also " + "be a comma-separated list."); + group.add_option("-X","--corrs",type="string", + help="correlation selection. Use comma-separated list of correlation indices."); + group.add_option("-S","--stations",type="string", + help="station (=antenna) selection. Use comma-separated list of station indices."), + group.add_option("-I","--ifrs",type="string", + help="interferometer selection. Use \"-I help\" to get help on selecting ifrs."); + group.add_option("-D","--ddid",type="string", + help="DATA_DESC_ID selection. Single number, or comma-separated list."); + group.add_option("-F","--field",type="string", + help="FIELD_ID selection. Single number, or comma-separated list."); + group.add_option("-Q","--taql",dest="taql",type="str", + help="additional TaQL selection to restrict subset."); + parser.add_option_group(group); + + group = OptionGroup(parser,"Selection by data value"); + group.add_option("--above",metavar="X",type="float", + help="select on abs(data)>X"); + group.add_option("--below",metavar="X",type="float", + help="select on abs(data) MS is %s"%msname; + print " %d antennas: %s"%(len(ants)," ".join(ants)); + print " %d DATA_DESC_ID(s): "%len(spwids); + for i,(spw,pol) in enumerate(zip(spwids,polids)): + print " %d: %.3f MHz, %d chans x %d correlations"%(i,ref_freq[spw]*1e-6,nchan[spw],len(corrs[pol,:])); + print " %d field(s): %s"%(len(fields),", ".join(["%d: %s"%ff for ff in enumerate(fields)])); + if not flagger.has_bitflags: + print "No BITFLAG/BITFLAG_ROW columns in this MS. Use the 'addbitflagcol' command to add them."; + else: + names = flagger.flagsets.names(); + if names: + print " %d flagset(s): "%len(names); + for name in names: + mask = flagger.flagsets.flagmask(name); + print " '%s': %d (0x%02X)"%(name,mask,mask); + else: + print " No flagsets."; + print ""; + if options.flag or options.unflag or options.fill_legacy or options.remove: + print "-l/--list was in effect, so all other options were ignored."; + sys.exit(0); + + # --flag/--unflag/--remove implies '-g all' by default, '-g -' skips the fill-legacy step + if options.flag or options.unflag or options.remove: + if options.fill_legacy is None: + options.fill_legacy = 'all'; + elif options.fill_legacy == '-': + options.fill_legacy = None; + + # if no other actions supplied, enable stats + if not (options.flag or options.unflag or options.fill_legacy): + statonly = True; + else: + statonly = False; + + # convert all the various FLAGS to flagmasks (or Nones) + for opt in 'data_flagmask','flagmask','flagmask_all','flagmask_none','flag','unflag','fill_legacy': + value = getattr(options,opt); + try: + flagmask = flagger.lookup_flagmask(value,create=(opt=='flag' and options.create)); + except Exception,exc: + msg = str(exc); + if opt == 'flag' and not options.create: + msg += "\nPerhaps you forgot the -c/--create option?"; + error(msg); + setattr(options,opt,flagmask); + + # clear the legacy flag itself from fill_legacy, otherwise it can have no effect + if options.fill_legacy is not None: + options.fill_legacy &= ~Flagger.LEGACY; + + # + # -r/--remove: remove flagsets + # + if options.remove is not None: + if options.flag or options.unflag: + error("Can't combine -r/--remove with --f/--flag or -u/--unflag."); + # get list of names and flagmask to remove + remove_names = options.remove.split(","); + names_found = []; + names_not_found = []; + flagmask = 0; + for name in remove_names: + try: + flagmask |= flagger.flagsets.flagmask(name); + names_found.append(name); + except: + names_not_found.append(name); + if not names_found: + if names_not_found: + error("No such flagset(s): %s"%",".join(names_not_found)); + error("No flagsets to remove"); + # unflag flagsets + if names_not_found: + print "===> WARNING: flagset(s) not found, ignoring: %s"%",".join(names_not_found); + print "===> removing flagset(s) %s"%",".join(names_found); + print "===> and clearing corresponding flagmask %s"%Flagger.flagmaskstr(flagmask); + if options.fill_legacy is not None: + print "===> and filling FLAG/FLAG_ROW using flagmask %s"%Flagger.flagmaskstr(options.fill_legacy); + flagger.xflag(unflag=flagmask,fill_legacy=options.fill_legacy); + flagger.flagsets.remove_flagset(*names_found); + sys.exit(0); + + # parse subset options + subset = parse_subset_options(options); + + # at this stage all remaining options are handled the same way + flagstr = unflagstr = legacystr = None; + if options.flag is not None: + flagstr = flagger.flagmaskstr(options.flag); + print "===> flagging with flagmask %s"%flagstr; + if options.unflag is not None: + unflagstr = flagger.flagmaskstr(options.unflag); + print "===> unflagging with flagmask %s"%unflagstr; + if options.fill_legacy is not None: + legacystr = flagger.flagmaskstr(options.fill_legacy); + print "===> filling legacy flags with flagmask %s"%legacystr; + + # do the job + totrows,sel_nrow,sel_nvis,nvis_A,nvis_B,nvis_C = \ + flagger.xflag(flag=options.flag,unflag=options.unflag,fill_legacy=options.fill_legacy,**subset); + + # print stats + if statonly: + print "===> No actions were performed. Showing the result of your selection:" + else: + print "===> Flagging stats:"; + rpc = 100.0/totrows; + print "===> MS size: %8d rows"%totrows; + print "===> Data selection: %8d rows, %8d visibilities (%.3g%% of MS rows)"%(sel_nrow,sel_nvis,sel_nrow*rpc); + if legacystr: + print "===> (over which legacy flags were filled using flagmask %s)"%legacystr; + + percent = 100.0/sel_nvis; + if options.channels or options.corrs: + print "===> Chan/corr slicing reduces this to %8d visibilities (%.3g%% of selection)"%(nvis_A,nvis_A*percent); + if not (options.flagmask is None and options.flagmask_all is None and options.flagmask_none is None): + print "===> Flag selection reduces this to %8d visibilities (%.3g%% of selection)"%(nvis_B,nvis_B*percent); + if options.above is not None or options.below is not None or options.fm_above is not None or options.fm_below is not None: + print "===> Data clipping reduces this to %8d visibilities (%.3g%% of selection)"%(nvis_C,nvis_C*percent); + if unflagstr: + print "===> (which were unflagged using flagmask %s)"%unflagstr; + if flagstr: + print "===> (which were flagged using flagmask %s)"%flagstr; + + flagger.close(); \ No newline at end of file diff --git a/merge-ms.py b/merge-ms.py new file mode 100755 index 0000000..7d4ee19 --- /dev/null +++ b/merge-ms.py @@ -0,0 +1,122 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import sys +import os.path +import re +import time + +time_start = time.time(); + +from Owlcat import table,tablecopy,tableexists,tabledelete +import Owlcat.Console; + +progress = Owlcat.Console.Reporter(timestamp=True); + +SKIP_COLUMNS = set(("FLAG_CATEGORY","NFRA_AVERAGEDATA","PEELED_DATA","UVMODEL_DATA")); + +if __name__ == "__main__": + + # setup some standard command-line option parsing + # + from optparse import OptionParser + parser = OptionParser(usage="""%prog: [options] MS_output MS_in1 MS_in2 ...""", + description="Concatenates several MSs into one. This is done on a row-by-row basis " + "with no sanity checking, so it's up to you to insure that the input MSs actually " + "belong together. All subtables and keywords of the output MS are inherited from " + "the first input MS, unless the --renumber-spws option is in effect, in which case " + "the SPECTRAL_WINDOW and DATA_DESCRIPTION tables are also merged (and spwids and " + "ddids renumbered.)"); + parser.add_option("-f","--force",dest="force",action="store_true", + help="overwrites output MS if it already exists"); + parser.add_option("-s","--renumber-spws",dest="renumber",action="store_true", + help="treat each MS as a separate spectral window"); + + (options,msnames) = parser.parse_args(); + + if len(msnames) < 3: + parser.error("Insufficient number of arguments. Use '-h' for help."); + + msout = msnames[0]; + msins = msnames[1:]; + + if tableexists(msout): + if not options.force: + print "Output MS",msout,"already exists. Use the -f switch to overwrite."; + sys.exit(1); + + if not options.force: + print "Will create merged MS %s from %d input MSs:"%(msout,len(msins)); + for ms in msins: + print " ",ms; + if options.renumber: + print "Each input MS will be put into a separate spectral window, spws will be renumbered."; + else: + print "Spectral windows will not be renumbered."; + if raw_input("Proceed (y/n)? ").strip().upper()[0] != "Y": + print "Aborted by user."; + sys.exit(1); + + if tableexists(msout) and options.force: + progress("Deleting existing copy of %s"%msout); + tabledelete(msout); + + # copy first MS to output as-is + progress("Copying %s to %s"%(msins[0],msout)); + tablecopy(msins[0],msout,deep=True); + + # open output for writing + tab0 = table(msout,readonly=False); + + # get the number of DDIDs and SPWIDs + if options.renumber: + ddid_tab0 = table(tab0.getkeyword("DATA_DESCRIPTION"),readonly=False); + spw_tab0 = table(tab0.getkeyword("SPECTRAL_WINDOW"),readonly=False); + + for msname in msins[1:]: + tab = table(msname); + nr0 = tab0.nrows(); + progress("Current row count: %d"%nr0); + progress("Merging in %d rows from %s"%(tab.nrows(),msname)); + tab0.addrows(tab.nrows()); + for col in tab.colnames(): + if col not in SKIP_COLUMNS: + data = tab.getcol(col); + # if renumbering DDIDs, increment the DDID column by the # of rows in the DDID table -- + # this ensures uniqueness of DDIDs + if options.renumber and col == "DATA_DESC_ID": + data += ddid_tab0.nrows(); + progress.overprint("Writing column %s, shape %s"%(col,data.shape)); + tab0.putcol(col,data,nr0); + # if renumbering, need to concatenate the DDID and SPW tables + if options.renumber: + progress.overprint("Updating DATA_DESCRIPTION subtable"); + # append content of the DATA_DESCRIPTION table, while renumbering spectral window IDs + ddid_tab = table(tab.getkeyword("DATA_DESCRIPTION"),readonly=False); + nr0 = ddid_tab0.nrows(); + ddid_tab0.addrows(ddid_tab.nrows()); + for col in ddid_tab.colnames(): + data = ddid_tab.getcol(col); + if col == "SPECTRAL_WINDOW_ID": + data += spw_tab0.nrows(); + progress.overprint("Writing column %s, shape %s"%(col,data.shape)); + ddid_tab0.putcol(col,data,nr0); + # append content of the SPECTRAL_WINDOW + progress.overprint("Updating SPECTRAL_WINDOW subtable"); + spw_tab = table(tab.getkeyword("SPECTRAL_WINDOW"),readonly=False); + nr0 = spw_tab0.nrows(); + spw_tab0.addrows(spw_tab.nrows()); + for col in spw_tab.colnames(): + data = spw_tab.getcol(col); + spw_tab0.putcol(col,data,nr0); + + + progress.overprint("Closing output MS %s\n"%msout); + nr0 = tab0.nrows(); + tab0.close(); + progress("Wrote %d rows to output MS"%nr0); + if options.renumber: + progress("Output MS contains %d spectral windows and %d DDIDs"% + (spw_tab0.nrows(),ddid_tab0.nrows())); + ddid_tab0.close(); + spw_tab0.close(); diff --git a/owlcat-logo.jpg b/owlcat-logo.jpg new file mode 100644 index 0000000..ce3cfd4 Binary files /dev/null and b/owlcat-logo.jpg differ diff --git a/plot-ms.py b/plot-ms.py new file mode 100755 index 0000000..f70c19f --- /dev/null +++ b/plot-ms.py @@ -0,0 +1,620 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import os.path +import math +import sys +import re +import warnings + +import numpy +import numpy.ma +import Owlcat +import matplotlib.pyplot as pyplot + +from Owlcat import Parsing + + +# +# NB: these need to be constructed on-the-fly based on polarization info from the MS. +# Otherwise we don't support circular without kludging +# + +# Each plotter func is a tuple of +# label,description,callable,flags_only +# The 'callable' transforms data into plottables +# If 'flags_only' is False, callable(visibilities) should return an array of real plottables +# the same shape as the visibilities masked_array. +# If 'flags_only' is False, callable(flags,axis) should return an array of flag densities +# along the specified axis. The returned array should be reduced along 'axis'. +# +Plotters = [ + ("I","Stokes I",lambda data:(abs(data[...,0])+abs(data[...,3]))/2,False), + ("Q","Stokes Q",lambda data:(abs(data[...,0])-abs(data[...,3]))/2,False), + ("U","Stokes U",lambda data:(data[...,1].real+data[...,2].real)/2,False), + ("V","Stokes V",lambda data:(data[...,1].imag-data[...,2].imag)/2,False), + ("flags_I","I/Q flag density", + lambda flags,meanaxis: + (flags[...,0]|flags[...,3]).mean(meanaxis),True ), + ("flags_all","2x2 flag density", + lambda flags,meanaxis: + (flags[...,0]|flags[...,1]|flags[...,2]|flags[...,3]).mean(meanaxis),True ), +]; + +for icorr,corr in enumerate(("XX","XY","YX","YY")): + Plotters += [ + (corr,corr+" amplitude",lambda data,i=icorr:abs(data[...,i]),False ), + (corr+"phase",corr+" phase",lambda data,i=icorr:numpy.ma.masked_array(numpy.angle(data[...,i]), + data[...,i].mask),False ), + (corr+"r",corr+" real",lambda data,i=icorr:data[...,i].real,False ), + (corr+"i",corr+" imag",lambda data,i=icorr:data[...,i].imag,False ), + ("flags_"+corr,corr+" flag density", + lambda flags,meanaxis,i=icorr: + flags[...,i].mean(meanaxis),True ), + ]; + + +if __name__ == "__main__": + + # setup some standard command-line option parsing + # + from optparse import OptionParser,OptionGroup + parser = OptionParser(usage="""%prog: [options] MS column:plot [column:plot ...]""", + description="""Makes plots of a column in the MS. Plots are specified as "column:plottype". If +column is not given, it defaults to the previous column, or CORRECTED_DATA. Use --list-plots to get +more information on available plots. If no plots are specified, CORRECTED_DATA:I is plotted by default. +"""); + + parser.add_option("--list-plots",action="store_true", + help="list available plot types and exit"); + parser.add_option("-L","--channels",dest="freqslice",type="string", + help="channel selection: single number or start:end[:step] to select channels start through end-1, " + "or start~end[:step] to select channels start through end, with an optional stepping."); + parser.add_option("-T","--timeslots",dest="timeslice",type="string", + help="timeslot selection by number, same format as channels."); + parser.add_option("-I","--ifrs",dest="ifrs",type="string", + help="subset of interferometers to plot. Use \"-I help\" to get help on selecting ifrs."); + parser.add_option("-D","--ddid",dest="ddid",type="string", + help="DATA_DESC_ID to plot. Default is first DDID found in MS. Use comma-separated list " + "(or 'all') to plot multiple DDIDs."); + parser.add_option("-F","--field",dest="field",type="int", + help="FIELD_ID to plot. Default is all (this may produce some strange plots for multiple-field MSs.)"); + parser.add_option("-Q","--taql",dest="taql",type="string", + help="additional TaQL selection to restrict data subset."); + parser.add_option("-f","--flagmask",metavar="FLAGS",dest="flagmask",type="string", + help="flagmask to apply to the BITFLAG column to get flagged data points" + "default is ALL, use 0 to ignore flags."); + + group = OptionGroup(parser,"Arranging your plots"); + group.add_option("--x-time",dest="xaxis",action="store_const",const=0, + help="use time for X axis and average over channels (default)"); + group.add_option("--x-freq",dest="xaxis",action="store_const",const=1, + help="use frequency for X axis and average over timeslots"); + + PLOT,DDID,IFR = group_choices = [ "plot","ddid","ifr" ]; + group.add_option("-P","--page",metavar="|".join(group_choices),action="append", + type="choice",choices=group_choices, + help="put different plots, DDIDs and/or interferometers on " + "separate pages. Can be given multiple times, this will determine " + "page order. Default is to put DDIDs and plots on different pages, and to " + "stack interferometers."); + group.add_option("-S","--stack",metavar="|".join(group_choices),action="append", + type="choice",choices=group_choices, + help="stack different plots, DDIDs and/or interferometers on " + "a single page. Can be given multiple times, this will determine " + "stacking order."); + choices = group_choices[1:]; + group.add_option("-A","--average",metavar="|".join(choices),action="append", + type="choice",choices=choices, + help="average DDIDs and/or interferometers together. Use twice " + "to average both."); + group.add_option("--ppp",dest="ppp",metavar="N",type="int", + help="maximum number of plots to stack per page. Default is 120."); + group.add_option("--offset",dest="offset",metavar="X",type="float", + help="vertical offset between stacked plots, in units of the median stddev. " + "Default is 10. NB: for flag-density plots, unit is 0.11*maxdensity, so the " + "standard offset of 10 produces plots spaced at 110% of the max value."); + parser.add_option_group(group); + + group = OptionGroup(parser,"Plot labels"); + group.add_option("--no-label-plot",dest="label_plot",action="store_false", + help="do not include plot type in labels (when stacking)"); + group.add_option("--no-label-ddid",dest="label_ddid",action="store_false", + help="do not include DDID in labels (when stacking)"); + group.add_option("--no-label-ifr",dest="label_ifr",action="store_false", + help="do not include IFR in labels (when stacking)"); + group.add_option("--no-label-baseline",dest="label_baseline",action="store_false", + help="do not include baseline length in labels"); + group.add_option("--no-label-mean",dest="label_mean",action="store_false", + help="do not include mean value in labels"); + group.add_option("--no-label-stddev",dest="label_stddev",action="store_false", + help="do not include standard deviation in labels"); + parser.add_option_group(group); + + group = OptionGroup(parser,"Output options"); + group.add_option("--title",dest="title",type="string",action="append", + help="provide a custom plot title. When you have multiple pages, you may " + "give this option serveral times."); + group.add_option("-o","--output",metavar="FILE",dest="output",type="string", + help="plot output to FILE. If not specified, plot will be " + "shown in a window. File format is determined by extension, see matplotlib documentation " + "for supported formats. At least .png, .pdf, .ps, .eps and .svg are supported."); + group.add_option("-r","--resolution",dest="resolution",type="int",metavar="DPI", + help="plot resolution for output to FILE (default is 100)"); + group.add_option("--size",metavar="WxH",dest="figsize",type="string", + help="set figure size, in cm. Default is '%default'."); + group.add_option("--papertype",dest="papertype",type="string", + help="set paper type (for .ps output only.) Prefix with '/' for " + "landscape mode. Default is '%default', but can also use e.g. 'letter', 'a3', etc."); + parser.add_option_group(group); + + + parser.set_defaults(output="",xaxis=0,ddid='first',field=None, + resolution=100,ppp=120,papertype='a4',figsize="21x29",offset=10, + flag_mask=None, + label_plot=True, + label_ddid=True, + label_ifr=True, + label_baseline=True, + label_mean=True, + label_stddev=True, + page=[],stack=[],average=[]); + + (options,args) = parser.parse_args(); + + # print help on plotters + if options.list_plots: + print "Available plot types:\n" + for p in Plotters: + print " %-16s%s"%(p[0],p[1]); + print """ +By default, specifying "DATA:XX" is the same as "DATA:XX.mean", producing a plot of mean |XX| +in time or frequency (depending on choice of X axis). Use "DATA.mean:XX" to plot |mean XX| +instead (i.e. mean visibilities, not mean amplitudes!) Use "DATA:XX.stddev" or "DATA.stddev:XX" +to plot standard deviations. Other options are ".sum", ".max" and ".min". + +Note that flag-density plots do not require a data column to be specified, since the +BITFLAG/FLAG columns are shared among all data columns. +"""; + sys.exit(0); + + # turn list of plotters into dict + Plotters = dict([(p[0],p[1:]) for p in Plotters]); + + # figure out plot arrangements, insert defaults for things that are not specified explicitly + allopt = options.page + options.stack + options.average; + if PLOT not in allopt: + options.page.insert(0,PLOT); + if DDID not in allopt: + options.page.insert(0,DDID); + if IFR not in allopt: + options.stack.insert(0,IFR); + # check that there's no double entries + for what in [PLOT,DDID,IFR]: + if len([opt for opt in allopt if opt == what ]) > 1: + parser.error("Conflicting --page/--stack/--average options for '%s'."); + + # get figsize + try: + w,h = map(int,options.figsize.split('x',1)); + figsize = (w*10,h*10); + except: + parser.error("Invalid --size setting '%s'"%options.figsize); + + # get MS name + if not args: + parser.error("MS not specified. Use '-h' for help."); + msname = args[0] + print "===> Attaching to MS %s"%msname; + ms = Owlcat.table(msname); + + # get flagmask, use a Flagger for this + import Owlcat.Flagger + from Owlcat.Flagger import Flagger + if options.flagmask == "0": + flagmask = 0; + print "===> Flagmask 0, ignoring all flags"; + elif options.flagmask is not None: + flagger = Flagger(msname); + flagmask = flagger.lookup_flagmask(options.flagmask); + flagger.close(); + flagger = None; + print "===> Flagmask is %s (you specified '%s')"%(Flagger.flagmaskstr(flagmask),options.flagmask); + if not flagmask&Flagger.LEGACY: + print "===> NB: legacy FLAG/FLAG_ROW columns will be ignored with this flagmask"; + else: + flagmask = Flagger.BITMASK_ALL|Flagger.LEGACY; + print "===> Using all flags"; + + # parse slice specs + try: + freqslice = Parsing.parse_slice(options.freqslice); + print "===> Channel selection is ",freqslice; + except: + parser.error("Invalid -L/--channels option. Start:end[:step] or start~end[:step] expected."); + try: + timeslice = Parsing.parse_slice(options.timeslice); + print "===> Timeslot selection is ",timeslice; + except: + parser.error("Invalid -T/--timeslots option. Start:end[:step] start~end[:step] expected."); + + # make list of plots + plots = []; + has_flagplots = False; + column0 = "CORRECTED_DATA"; + if column0 not in ms.colnames(): + column0 = "DATA"; + # go through list of arguments, or default list + for arg in (args[1:] or ["I"]): + # parse as "[column[.reduce]:]plottype[.reduce]" + m = re.match('^(\w+)(\.(\w+))?(:(\w+)(.(\w+))?)?$',arg); + if not m: + parser.error("'%s': invalid argument"%arg); + g = m.groups(); + # groups 0-2 are the first field, 3-6 are the second field + # if 3 didn't match, then we only have one field, which we treat as a plottype + if g[3]: + column,datareduce,plot,plotreduce = g[0],g[2],g[4],g[6]; + else: + column,datareduce,plot,plotreduce = None,None,g[0],g[2]; + # the reduction function applies to visibilities or to plottable values + if plotreduce and datareduce: + parser.error("'%s': can't specify two .funcs"%arg); + # check that function is valid + for func in datareduce,plotreduce: + if func not in [None,'mean','std','min','max','sum','product' ]: + parser.error("'%s': '%s' is not a valid reduction function"%(arg,func)); + # lookup plotter + plotdesc,func,is_flagplot = Plotters.get(plot,(None,None,None)); + if func is None: + parser.error("'%s': unknown plot type '%s'"%(arg,plot)); + # more sanity checks + if is_flagplot: + if plotreduce or datareduce: + parser.error("'%s': can't use '.%s' with flagplots"%(arg,(plotreduce or datareduce))); + if column: + parser.error("'%s': can't specify a column for flagplots"%arg); + else: + column0 = column or column0; + if column0 not in ms.colnames(): + parser.error("MS does not contain a %s column"%column); + # set default reduction, if nothing else is specified + if not plotreduce and not datareduce: + plotreduce = 'mean'; + plotdesc = " ".join(filter(bool,[datareduce,column0,plotreduce,plotdesc])); + # add to list of plots + has_flagplots |= is_flagplot; + plots.append((plot,plotdesc,func,is_flagplot,column0,datareduce,plotreduce)); + + # collect applicable TaQL queries here + taqls = []; + + # get IFR set + import Meow.IfrSet + ifrset = Meow.IfrSet.from_ms(ms); + tot_ifrs = len(ifrset.ifrs()); + if options.ifrs: + if options.ifrs == "help": + # print help string, but trim away RTF tags + print re.sub("<[^>]+>","",ifrset.subset_doc).replace("<","<").replace(">",">"); + sys.exit(0); + ifrset = ifrset.subset(options.ifrs); + taqls.append(ifrset.taql_string()); + print "===> Selected %d of %d interferometers "%(len(ifrset.ifrs()),tot_ifrs); + else: + print "===> Selected all %d interferometers "%tot_ifrs; + + # select DDIDs + ddid_tab = Owlcat.table(ms.getkeyword('DATA_DESCRIPTION')); + # default ('first') is to use first DDID in MS + ddid_str = options.ddid.strip().lower(); + if ddid_str == 'first': + ddids = [ ms.getcol('DATA_DESC_ID',0,1)[0] ]; + print "===> Using first DATA_DESC_ID (%d)"%ddids[0]; + # else see if it's a single int + elif re.match('^\d+$',ddid_str): + ddids = [ int(options.ddid) ]; + # else parse as list of ints, or 'all'. In this case we need extra info from the tables + else: + if ddid_str == "all": + ddids = range(ddid_tab.nrows()); + else: + try: + ddids = map(int,ddid_str.split(',')); + except: + parser.error("Invalid -D/--ddid option: %s"%ddid_str); + + # Get data shapes by reading subtables. I used to just use: + # datashape = [ ms.nrows() ] + list(ms.getcoldesc('DATA')['shape']); + # but this is not robust, since some MSs will not have a fixed-shape DATA column. So, read the subtables. + spwids = ddid_tab.getcol('SPECTRAL_WINDOW_ID'); + polids = ddid_tab.getcol('POLARIZATION_ID'); + corrs = Owlcat.table(ms.getkeyword('POLARIZATION')).getcol('CORR_TYPE'); + spw_tab = Owlcat.table(ms.getkeyword('SPECTRAL_WINDOW')); + ref_freq = spw_tab.getcol('REF_FREQUENCY'); + nchan = spw_tab.getcol('NUM_CHAN'); + + # select field + if options.field is not None: + taqls.append("FIELD_ID==%d"%options.field); + print "===> Selecting FIELD_ID %d"%options.field; + + # select TaQL + if options.taql: + taqls.append(options.taql); + print "===> Additional TaQL query is \"%s\""%options.taql; + + # apply accumulated selection + if taqls: + ms = ms.query("( " + " ) && ( ".join(taqls) + " )"); + print "===> Selected %d rows from MS"%ms.nrows(); + if not ms.nrows(): + print """MS selection is empty. You may have specified it incorrectly: please check your +DATA_DESC_ID (option -D/--ddid), field (-F/--field), interferometer subset (-I/--ifrs) +and/or TaQL query (-Q/--taql) options. Or was your MS empty to begin with?"""; + sys.exit(1); + + # get axis indices + xaxis = options.xaxis; + meanaxis = 1 if xaxis==0 else 0; + + # figure out label format + labels = []; + if options.label_plot and PLOT in options.stack: + labels.append("%(plot)s"); + if options.label_ddid and DDID in options.stack: + labels.append("d#%(ddid)d"); + if options.label_ifr and IFR in options.stack: + labels.append("%(ifr)s"); + if options.label_baseline and IFR in options.stack: + labels.append("%(baseline)dm"); + if options.label_mean: + labels.append("mean=%(mean).3g"); + if options.label_stddev: + labels.append("std=%(stddev).3g"); + label_format = " ".join(labels); + + # Make list of trackkeys + # A trackkey is (nplot,ddid,ifr); ddid or ifr is None if averaging + average_ddids = DDID in options.average; + average_ifrs = IFR in options.average; + # ranges for each trackkey + keyranges = [ range(len(plots)), + [None] if average_ddids else ddids, + [None] if average_ifrs else [ (px[0],qx[0]) for px,qx in ifrset.ifr_index() ] ]; + + import Owlcat.Plotting + + # A PlotCollection holds one set of plot tracks. + # This is a dict of trackkey:PC. + # PC's are shared by multiple tracks according to the page/stack/average settings. + plotcolls = {}; + + # Helper function to get a PC associated with the given trackkey + # Creates one if it doesn't exist, and inserts it into the plotcolls dict + # for every other track in the stacked set. + def get_plot_collection (track): + pc = plotcolls.get(track); + if pc: + return pc; + # make new PlotCollection + pc = plotcolls[track] = Owlcat.Plotting.PlotCollection(options); +# print "New pc for",track; + # associate it with every track that will go on it, by looping over stacked keys + for plot1 in (keyranges[0] if PLOT in options.stack else [track[0]]): + for ddid1 in (keyranges[1] if DDID in options.stack else [track[1]]): + for ifr1 in (keyranges[2] if IFR in options.stack else [track[2]]): + plotcolls[plot1,ddid1,ifr1] = pc; +# print "also for ",(plot1,ddid1,ifr1); + return pc; + + # set of IFRs for which (non-flagged) data is actually found + active_ifrs = set(); + labelattrs = {}; # dict of plot label attrs, updated in each loop + + # now loop over DDIDs + for ddid in ddids: + labelattrs['ddid'] = ddid; + + subms = ms.query("DATA_DESC_ID==%d"%ddid); + datashape = [ subms.nrows(),nchan[spwids[ddid]],len(corrs[polids[ddid]]) ]; + + print "===> Processing DATA_DESC_ID %d (%d MHz): %dx%d by %d rows"%( + ddid,round(ref_freq[spwids[ddid]]*1e-6),datashape[1],datashape[2],datashape[0]); + + if not subms.nrows(): + continue; + # get boolean flags + if flagmask&Flagger.LEGACY: + flagcol = subms.getcol('FLAG'); + # merge in FLAG_ROW column + flagrowcol = subms.getcol('FLAG_ROW'); + flagcol |= flagrowcol[:,numpy.newaxis,numpy.newaxis]; + else: + flagcol = numpy.zeros(datashape,bool); + # get bitflags + bitflags = flagmask&Flagger.BITMASK_ALL; + if bitflags: + if 'BITFLAG' in ms.colnames(): + bf = subms.getcol('BITFLAG'); + flagcol |= ((bf&bitflags)!=0); + if 'BITFLAG_ROW' in ms.colnames(): + bfr = subms.getcol('BITFLAG_ROW'); + flagcol |= ((bfr&bitflags)!=0)[:,numpy.newaxis,numpy.newaxis]; + flagcol = flagcol[:,freqslice,:]; + nf = flagcol.sum(); + print "===> %d of %d (%.2g%%) visibilities are flagged "%(nf,flagcol.size,(nf/float(flagcol.size))*100); + + # get antenna indices, and make dict per-ifr masks + a1,a2 = subms.getcol('ANTENNA1'),subms.getcol('ANTENNA2'); + ifr_rows = dict([ ((p,q),numpy.where((a1==p)&(a2==q))[0]) for (p,plab),(q,qlab) in ifrset.ifr_index() ]); + # visibility columns will be read into here, as demanded by the individual plots + datacols = {}; + + colname0 = None; + for iplot,(plottype,plotdesc,plotfunc,is_flagplot,colname,datareduce,plotreduce) in enumerate(plots): + if colname0 and colname and colname != colname0: + labelattrs['plot'] = "%s:%s"%(colname,plottype); + else: + labelattrs['plot'] = plottype; + colname0 = colname; + print "===> Making plot",plotdesc; + # read in data column, if not already read + if colname: + # do we already have a reduced version of the column cached? + dc = datacols.get((colname,datareduce),None); + if dc is None: + # no, do we have an unreduced version cached? + dc = datacols.get((colname,None),None); + if dc is None: + # ok, read & cache column + dc = subms.getcol(colname)[:,freqslice,:]; + dc = datacols[colname,None] = numpy.ma.masked_array(dc,flagcol,fill_value=complex(0,0)); + # reduce and cache reduced version + if datareduce: + dc = datacols[colname,datareduce] = getattr(dc,datareduce)(meanaxis); + + # split up into per-baseline arrays, and accumulate data tracks + for (p,plab),(q,qlab) in ifrset.ifr_index(): + idx = ifr_rows[p,q][timeslice]; + if not len(idx): + continue; + # this is the plot track ID. Things that we average over are replaced by None + track = (iplot,None if average_ddids else ddid,None if average_ifrs else (p,q)); + # update label attributes + labelattrs['baseline'] = round(ifrset.baseline(p,q)); + labelattrs['ifr'] = ifrset.ifr_label(p,q); + # + if is_flagplot: + d1 = numpy.ma.masked_array(plotfunc(flagcol[idx,...],meanaxis)); + else: + d1 = plotfunc(dc[idx,...]); + if plotreduce: + d1 = getattr(d1,plotreduce)(meanaxis); + d1.fill_value = 0; + # check that this data track is not fully flagged + if d1.mask.all(): + continue; + active_ifrs.add((p,q)); + # get a plot collection object, or make a new one + plotcoll = get_plot_collection(track); + # for flag-density plots, set the plot offset + if is_flagplot: + plotcoll.offset = max(plotcoll.offset,d1.max()*0.11*options.offset); + # add or update track + count = plotcoll.get_track_count(track); + # count !=0 indicates track is being updated (because we're averaging DDIDs or IFRs) + if count: + d0 = plotcoll.get_track_data(track); + if d0.shape != d1.shape: + warnings.warn("Shape mismatch between averaged DDIDs or IFRs, will ignore mismatched data."); + continue; + # add data to track. For means, update as (old*n + new)/(n+1), but if stddev + # reduction was in effect, update as sqrt((old**2)*n + new)/(n-1) + if 'std' in (plotreduce,datareduce): + d1 = numpy.sqrt(((d0**2)*count+d1**2)/(count+1)); + else: + d1 = (d0*count+d1)/(count+1); + # insert track data + labelattrs['mean'] = mean = d1.mean(); + labelattrs['stddev'] = std = d1.std(); + plotcoll.add_track(track,d1,count=count+1,mean=mean,stddev=std,label=label_format%labelattrs); + + # deallocate datacubes + datacols = flagcol = dc = None; + + # make list of active IFRS (as p,q pairs), sorted by baseline length + print "===> Found data for %d interferometers"%len(active_ifrs); + if not average_ifrs: + keyranges[2] = sorted(active_ifrs,lambda a,b:cmp((round(ifrset.baseline(*a)),a),(round(ifrset.baseline(*b)),b))); + + if not active_ifrs: + print "===> Nothing to be plotted. Check your data selection."; + sys.exit(0); + + + # Now, work out the order in which we plot the tracks and PlotCollections. + # this is determined by the options.page and options.stack order + + # Form up list of all possible track keys, in the order in which they are + # to be plotted. Note that after this first loop, the elements of each key are + # scrambled, because of the way we build up the list... + track_order = [[]]; + elorder = {}; + keyels = { PLOT:0,DDID:1,IFR:2 }; + for i,what in enumerate(options.average + options.stack + options.page): + elorder[what] = i; + track_order = [ to+[key] for key in keyranges[keyels[what]] for to in track_order ]; + # ...so unscamble them into proper order and turn them into tuples + track_order = [ ( key[elorder[PLOT]],key[elorder[DDID]],key[elorder[IFR]] ) for key in track_order ]; + + # figure out paper size + landscape = options.papertype[0] == "/"; + papertype = options.papertype[1:] if landscape else options.papertype; + + # now, split this up into PlotCollections and the tracks associated with each + pc_tracks = []; + pc0 = None; + for key in track_order: + pc = plotcolls[key]; + # start new collection + if pc is not pc0: + pc_tracks.append((pc,[key])); + pc0 = pc; + # else key belongs to previous collection + else: + pc_tracks[-1][1].append(key); + + # Do we have multiple plots, or multiple pages per plot? Start a page counter then. + if len(pc_tracks) > 1 or len(pc_tracks[0][1]) > options.ppp: + ipage = 0; + else: + ipage = None; + + # loop over plotcollectors + for iplotcoll,(plotcoll,keylist) in enumerate(pc_tracks): + # do we have a title for this plot on the command line? + if options.title and len(options.title) > iplotcoll: + title0 = options.title[iplotcoll]; + # no, make default title + else: + nplot,ddid,ifr = keylist[0]; + chans = freqslice.indices(datashape[1]); + if chans[1] == chans[0]+1: + chandesc = str(chans[0]); + elif chans[2] == 1: + chandesc = "%d~%d"%(chans[0],chans[1]); + else: + chandesc = "%d~%d step %d"%(chans[0],chans[1]-1,chans[2]); + title0 = msname; + if average_ddids: + title0 += " ddid %s"%"+".join(map(str,ddids)); + else: + title0 += " ddid %d"%ddid; + if options.freqslice: + title0 += " (chan %s)"%chandesc; + if average_ifrs: + title0 += " mean of %d ifrs"%len(active_ifrs); + elif IFR in options.page: + title0 += " ifr %s (%dm)"%(ifrset.ifr_label(*ifr),round(ifrset.baseline(*ifr))); + # plot description + if PLOT in options.page: + title0 += " "+plots[nplot][1]; + # split keylist if needed + for nkey0 in range(0,len(keylist),options.ppp): + nkey1 = min(nkey0+options.ppp,len(keylist)); + keys = keylist[nkey0:nkey1]; + # make filename, if saving file + savefile = options.output; + if ipage is not None: + if savefile: + basename,ext = os.path.splitext(savefile); + savefile = "%s.%d%s"%(basename,ipage,ext); + ipage += 1; + dual = (len(keys) > options.ppp/2); + plotcoll.make_figure(keys,save=savefile,suptitle=title0,figsize=figsize,offset_std=options.offset, + papertype=papertype,dpi=options.resolution,dual=dual,landscape=landscape); + + if not options.output: + from pylab import plt + plt.show(); diff --git a/restore-sky-model.py b/restore-sky-model.py new file mode 100755 index 0000000..1dc18a8 --- /dev/null +++ b/restore-sky-model.py @@ -0,0 +1,88 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import sys +import pyfits +import re +import os.path +import pyfits + + +if __name__ == '__main__': + # setup some standard command-line option parsing + # + from optparse import OptionParser + parser = OptionParser(usage="""%prog: [options] input_image newstar_sky_model [output_image]"""); + parser.add_option("-n","--num-sources",dest="nsrc",type="int",action="store", + help="Only restore the NSRC brightest sources"); + parser.add_option("-s","--scale",dest="fluxscale",metavar="FLUXSCALE[,N]",action="store", + help="rescale model fluxes by given factor. If N is given, rescale N brightest only."); + parser.add_option("-b","--beamsize",dest="beamsize",type="float",action="store", + help="restoring beam size (0 for a delta-function.)"); + parser.add_option("-p","--psf",dest="psf",action="store", + help="name of PSF file to be fitted, if beam size is not specified (default psf.fits)"); + parser.add_option("-f",dest="force",action="store_true", + help="overwrite output image even if it already exists"); + parser.add_option("-v","--verbose",dest="verbose",type="int",action="store", + help="set verbosity level (0 is silent, higher numbers mean more messages)"); + parser.set_defaults(n=0,fluxscale='1',psf='psf.fits',beam=-1); + + (options,rem_args) = parser.parse_args(); + + # get filenames + if len(rem_args) == 2: + input_image,skymodel = rem_args; + name,ext = os.path.splitext(input_image) + output_image = name+".restored"+ext; + elif len(rem_args) == 3: + input_image,skymodel,output_image = rem_args; + else: + parser.error("Insufficient number of arguments. Use -h for help."); + + # check for overwritten output + if os.path.exists(output_image) and not options.force: + print "File %s already exists, use the -f option to overwrite."%output_image; + sys.exit(1); + + from Timba.Contrib.OMS.SkyPuss.Tools import Imaging + from Timba.Contrib.OMS.SkyPuss import Import + + Imaging._verbosity.set_verbose(options.verbose); + + # read model and sort by apparent brightness + model = Import.importNEWSTAR(skymodel); + Imaging.dprintf(1,"Read %d sources from %s\n",len(model.sources),skymodel); + sources = sorted(model.sources,lambda a,b:cmp(b.brightness(),a.brightness())); + + # apply counts and flux scales + if options.nsrc: + sources = sources[:options.nsrc]; + Imaging.dprintf(1,"Using %d brightest sources\n",len(sources)); + + if options.fluxscale != '1': + if "," in options.fluxscale: + scale,n = options.fluxscale.split(",",1); + scale = float(scale); + n = int(n); + Imaging.dprintf(1,"Flux of %d brightest sources will be scaled by %f\n",n,scale); + else: + scale = float(options.fluxscale); + n = len(sources); + Imaging.dprintf(1,"Flux of all model sources will be scaled by %f\n",n,scale); + for src in sources[:n]: + src.flux.rescale(0.01); + + if options.beamsize >= 0: + gx = gy = options.beam; + grot = 0; + elif options.psf: + # fit the psf + gx,gy,grot = Imaging.fitPsf(options.psf); + + # read, restore, write + input_hdu = pyfits.open(input_image)[0]; + Imaging.restoreSources(input_hdu,sources,gx,gy,grot); + + Imaging.dprintf(1,"Writing output file %s\n",output_image); + input_hdu.writeto(output_image,clobber=True); + diff --git a/run-imager.sh b/run-imager.sh new file mode 100755 index 0000000..b265f9f --- /dev/null +++ b/run-imager.sh @@ -0,0 +1,176 @@ +# -*- coding: utf-8 -*- +# if config doesn't exist, strange things happen! + +# if dirty image is not made first, clean fails! +# Same message as when ddid_index is not set properly: +# Maximum of approximate PSF for field 1 = 0 : renormalizing to unity +# But not sure what state persists between making a dirty image and a clean one + +CONFFILE=imager.conf + +# Unset all variables beginning with img_ + + +# Default settings +img_lwimager=lwimager +img_data=DATA +img_ifrs="*" +img_weight=natural +img_spwid=0 +img_field=0 +img_size=512/60 +img_mode=channel +img_stokes=I +img_padding=1.2 +img_flux_scale=2 + +img_channels="" +img_img_channels="" + +img_remove_img=1 + +img_name_dirty=dirty +img_name_restored=restored +img_name_model=model +img_name_residual=residual + +img_cachesize=512 + +img_oper=image + +img_niter=1000 +img_gain=.1 +img_threshold=0Jy + +# load conf file +if [ -f $CONFFILE ]; then + source $CONFFILE +fi + +trial=true +confirm=true + +# overwrite with command line +while [ "$1" != "" ]; do + if [ "${1%=*}" != "$1" ]; then + eval "img_$1" + elif [ "$1" == "-nr" ]; then + img_remove_img=0 + elif [ "$1" == "-trial" ]; then + trial="echo 'Trial mode only, exiting'; exit 1"; + elif [ "$1" == "-i" ]; then + confirm="echo 'Press Enter to continue with these settings, or Ctrl+C to stop'; read"; + fi + shift +done + +# expand vars in filenames +img_name_dirty=`eval echo $img_name_dirty` +img_name_restored=`eval echo $img_name_restored` +img_name_model=`eval echo $img_name_model` +img_name_residual=`eval echo $img_name_residual` + +# MS must be set +if [ -z "$img_ms" ]; then + echo "MS must be set with MS=name, or in $CONFFILE"; +fi + +# print final var settings +set | grep ^img_ | cut -c 5- + +# expand size into arcmin and npix +if [ "$img_arcmin" == "" -o "$img_npix" == "" ]; then + img_arcmin=${img_size#*/} + img_npix=${img_size%/*} +fi + +# print ifr settings +if [ "$img_ifrs" == "" ]; then + echo "Using all interferometers" +else + ifr_select="`ifr-subset-to-taql.sh $img_ms ${img_ifrs// /,}`" + echo "NB: IFR selection '$img_ifrs' applied, resulting in ${ifr_select%//*} ifrs" + ifr_select="${ifr_select#*//}" + if [ "$img_select" != "" ]; then + img_select="($img_select)&&($ifr_select)" + else + img_select="$ifr_select" + fi + echo "select=$img_select" +fi + +eval $confirm + +CELLSIZE=`python -c "print $img_arcmin*60./float($img_npix)"` + +# make dirty image +make_image () +{ + # collect mandatory arguments + cmd="$img_lwimager ms=$img_ms data=$img_data operation=$img_oper + stokes=$img_stokes mode=$img_mode weight=$img_weight + npix=$img_npix cellsize=${CELLSIZE}arcsec + spwid=$img_spwid field=$img_field padding=$img_padding cachesize=$img_cachesize + " + + # add optional arguments + if [ "$img_weight" == "radial" -a "$img_taper" != "" ]; then + cmd="$cmd filter=${img_taper}arcsec,${img_taper}arcsec,0.000000deg" + fi + if [ "$img_channels" != "" ]; then + declare -a chans=(${img_channels//,/ }) + cmd="$cmd chanmode=channel nchan=${chans[0]} chanstart=${chans[1]} chanstep=${chans[2]}"; + fi + if [ "$img_img_channels" != "" ]; then + declare -a chans=(${img_img_channels//,/ }) + cmd="$cmd img_nchan=${chans[0]} img_chanstart=${chans[1]} img_chanstep=${chans[2]}"; + fi + if [ "$img_select" != "" ]; then + cmd="$cmd select='$img_select'" + fi + + # setup output files + if [ "$img_oper" == "image" ]; then + imgname="$img_name_dirty.img" + imgname_fits="$img_name_dirty.fits" + cmd="$cmd image=$imgname" + rm -fr $imgname_fits $imgname &>/dev/null + echo "Making dirty image: " $cmd + else + restored="$img_name_restored.img"; + model="$img_name_model.img"; + residual="$img_name_residual.img"; + restored_fits="$img_name_restored.fits"; + model_fits="$img_name_model.fits"; + residual_fits="$img_name_residual.fits"; + rm -fr $restored_fits $model_fits $residual_fits $restored $model $residual &>/dev/null + cmd="$cmd model=$model restored=$restored residual=$residual + niter=$img_niter gain=$img_gain threshold=$img_threshold + " + echo "Making clean image: " $cmd + fi + + eval $trial + if eval time $cmd; then + if [ "$img_remove_img" != "0" ]; then + echo "Success, all *.img files will be removed after conversion to FITS" + remove=rm + else + echo "Success" + remove=true + fi + if [ "$img_oper" == "image" ]; then + image2fits in=$img_flux_scale\*\\$imgname out=$imgname_fits && $remove -fr $imgname + else + image2fits in=$img_flux_scale\*\\$model out=$model_fits && $remove -fr $model + image2fits in=$img_flux_scale\*\\$residual out=$residual_fits && $remove -fr $residual + image2fits in=$img_flux_scale\*\\$restored out=$restored_fits && $remove -fr $restored + fi + return 0 + else + echo "Imager failed"; + return 1 + fi +} + +make_image diff --git a/split-ms-spw.py b/split-ms-spw.py new file mode 100755 index 0000000..f99c704 --- /dev/null +++ b/split-ms-spw.py @@ -0,0 +1,65 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import sys +import os.path +import re +import time + +from Owlcat import table,tablecopy,tableexists,tabledelete + +if __name__ == "__main__": + + # setup some standard command-line option parsing + # + from optparse import OptionParser + parser = OptionParser(usage="""%prog: [options] MS""", + description="Splits MS into one or more output MSs by DATA_DESC_ID (which typically corresponds" + "to spectral window.)"); + parser.add_option("-D","--ddid",type="int", + help="specific DATA_DESC_ID to extract. If not given, then all DDIDs are extracted, each into a separate MS."); + parser.add_option("-o","--output",type="string", + help="Name of output MS. Use '%(ms)s' to insert basename of input MS. Use '%(ddid)d' to insert DDID number. Use '%(msext)s' to insert extension of input MS. Default is %default"); + parser.add_option("-f","--force",dest="force",action="store_true", + help="overwrites output MS if it already exists"); + + parser.set_defaults(ddid=-1,output="%(ms)s_spw%(ddid)d%(msext)s"); + + (options,msname) = parser.parse_args(); + + if len(msname) != 1 : + parser.error("Incorrect number of arguments. Use '-h' for help."); + + # open input MS + msname = msname[0]; + print "Reading input MS %s"%msname; + if not tableexists(msname): + parser.error("MS %s not found."%msname); + ms = table(msname); + # check DDID option + num_ddids = table(ms.getkeyword("DATA_DESCRIPTION")).nrows(); + print "MS contains %d DATA_DESC_IDs"%num_ddids; + if options.ddid < 0: + ddids = range(num_ddids); + elif options.ddid >= num_ddids: + parser.error("DATA_DESC_ID %d is out of range"%options.ddid); + else: + ddids = [ options.ddid ]; + + # setup outputs + msname,msext = os.path.splitext(os.path.basename(os.path.normpath(msname))); + + for ddid in ddids: + msout = options.output%dict(ms=msname,msext=msext,ddid=ddid); + print "Extracting DATA_DESC_ID %d into MS %s"%(ddid,msout); + # delete if exists + if tableexists(msout): + if not options.force: + parser.error("Output MS %s already exists. Use the -f switch to overwrite."%msout); + print "Deleting existing copy of %s"%msout; + tabledelete(msout); + # extract + ms1 = ms.query("DATA_DESC_ID==%d"%ddid); + ms1.copy(msout,deep=True); + ms1.close(); + diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/.purrlock b/tutorial/Owlcat-plotms-tutorial.purrlog/.purrlock new file mode 100755 index 0000000..604437a --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/.purrlock @@ -0,0 +1 @@ +jakob:10687 \ No newline at end of file diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/dirconfig b/tutorial/Owlcat-plotms-tutorial.purrlog/dirconfig new file mode 100644 index 0000000..67ad61c --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/dirconfig @@ -0,0 +1,3 @@ +[~/Waterhole/contrib/OMS/Owlcat/tutorial] +watching = 1 + diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/index.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/index.html new file mode 100644 index 0000000..1035450 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/index.html @@ -0,0 +1,23 @@ + + Introducing Owlcat

Introducing Owlcat

+

Logged on 02/27/2010 02:12:41 AM

+ + + +

OWLCAT stands for Oleg's W--- Lightweight Calibration & Analysis Toolkit. The W can represent Wisely, Wantonly, Woefully, Weirdly, Wonderfully, Wickedly, Wobbly, Wankingly, Warp 9, or Whatever suits your mood. Owlcat is a bunch of scripts I wrote to help me reduce data. Of these, plot-ms.py, flag-ms.py and plot-parms.py are probably the most useful.

+

Owlcat is (currently) in my Waterhole, so if you want to use it, you must first check it out:

+

$ svn co svn://lofar9.astron.nl/var/svn/repos/trunk/Waterhole/contrib/OMS/Owlcat

+

This will create a directory called Owlcat. Now, make sure this directory is both in your PYTHONPATH and PATH.

+

Before you go any further, please make sure you've got Pylab/Matplotlib installed (Ubuntu package python-matplotlib), and have a copy of the MeqTrees Cattery scripts (svn co svn://lofar9.astron.nl/var/svn/repos/trunk/Frameworks/Cattery. No working MeqTrees required, just a copy of the Cattery python scripts.) Owlcat will not work without these.

+ + +

Data products

+ + + + + + + + +
owlcat-logo.jpg

owlcat-logo.jpg: The Owlcat logo

diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/index.include.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/index.include.html new file mode 100644 index 0000000..e1eac1c --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/index.include.html @@ -0,0 +1,22 @@ + +
+

Introducing Owlcat

+

Logged on 02/27/2010 02:12:41 AM

+ + + +

OWLCAT stands for Oleg's W--- Lightweight Calibration & Analysis Toolkit. The W can represent Wisely, Wantonly, Woefully, Weirdly, Wonderfully, Wickedly, Wobbly, Wankingly, Warp 9, or Whatever suits your mood. Owlcat is a bunch of scripts I wrote to help me reduce data. Of these, plot-ms.py, flag-ms.py and plot-parms.py are probably the most useful.

+

Owlcat is (currently) in my Waterhole, so if you want to use it, you must first check it out:

+

$ svn co svn://lofar9.astron.nl/var/svn/repos/trunk/Waterhole/contrib/OMS/Owlcat

+

This will create a directory called Owlcat. Now, make sure this directory is both in your PYTHONPATH and PATH.

+

Before you go any further, please make sure you've got Pylab/Matplotlib installed (Ubuntu package python-matplotlib), and have a copy of the MeqTrees Cattery scripts (svn co svn://lofar9.astron.nl/var/svn/repos/trunk/Frameworks/Cattery. No working MeqTrees required, just a copy of the Cattery python scripts.) Owlcat will not work without these.

+ + +

Data products

+ + + + + + +
owlcat-logo.jpg

owlcat-logo.jpg: The Owlcat logo

\ No newline at end of file diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/owlcat-logo.jpg b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/owlcat-logo.jpg new file mode 100644 index 0000000..ce3cfd4 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/owlcat-logo.jpg differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/owlcat-logo.jpg.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/owlcat-logo.jpg.purr-products/thumb.png new file mode 100644 index 0000000..c919794 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-021241/owlcat-logo.jpg.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/index.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/index.html new file mode 100644 index 0000000..ff10555 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/index.html @@ -0,0 +1,103 @@ + + Different kinds of plots

Different kinds of plots

+

Logged on 02/27/2010 02:23:58 AM

+ + + +

Currently, the most versatile tool in the toolkit is the plot-ms.py script. To try it out, arm yourself with an MS. I will provide some examples based on the 3C147_spw0.MS, which you can download here:

+

http://www.astron.nl/meqwiki-data/users/oms/3C147-Calibration-Tutorial/3C147_spw0.MS.tgz

+

For starters, try

+

$ plot-ms.py -h

+

This just to let you know that there's many options, and I've tried to document them all. Don't bother reading them in detail, just type

+

$ plot-ms.py 3C147_spw0.MS

+

to make your first plot. This is a plot of the mean Stokes I (|XX|+|YY|) in the CORRECTED_DATA column, as a function of time, averaged over all frequencies. Note how data from individual baselines is stacked up and sorted by baseline length.

+

The plot is displayed in a window. Close the window to exit the script. To save the plot to a file, use the -o option:

+

$ plot-ms.py 3C147_spw0.MS -o plot1.png

+

Besides .png, you can also use .ps, .eps, .svg, and whatever else is supported by matplotlib -- see matplotlib documentation. For PostScript output, you can set paper size, plot size, etc -- see the help for details.

+

Now, this was a plot of mean Stokes I in CORRECTED_DATA. This is the default plot type made when nothing else is specified. To make different plots, you must specify them on the command line after the MS name, as "column:plottype". For example

+

$ plot-ms.py 3C147_spw0.MS DATA:XX DATA:YY

+

makes two plots, of mean |XX| and |YY| in the DATA column. Note the label at the top of the plot, which gives a bunch of details.

+

Instead of the mean |XX| in frequency, you can also plot standard deviations:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX.std

+

Other functions you can try are .min, .max and .sum. The default plot (plot 1) is equivalent to CORRECTED_DATA:I.mean. Finally, try:

+

$ plot-ms.py 3C147_spw0.MS CORRECTED_DATA.mean:I

+

and see if you can spot the difference. In plot 1, we plotted mean amplitudes, here we plot the amplitude of the (complex) mean visibility. This is an important distinction (for starters, the mean complex visibility always tends to be smaller, as it is a vector sum that tends to "wash out" due to the phase component.)! The positioning of ".mean" determines whether it applies to the column (i.e. visibilities), or to the amplitudes.

+

The previous plots were a function of time & averaged in frequency. To plot things as a function of frequency (averaged in time), use the --x-freq option. This will show you the time-frequency amplitude and phase bandpasses:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq

+

The gap in the plot is due to a flagged channel. To ignore flags (i.e. plot flagged data anyway), use the -F option. Note that this may impact the scale of your plot (which is chosen automatically based on the variation in the plotted data.)

+

Speaking of flags, another interesting plot type is the flag density plot. This shows the fraction of visibilities flagged per time slot or frequency channel. Use

+

$ plot-ms.py 3C147_spw0.MS flags_XX

+

to get a flag-density plot for XX (as a function of time.) Add the --x-freq option to plot it as a function of frequency. If you're flagging RFI, these kinds of plots are a nice way to see the bad timeslots and frequency channels.

+

Finally, note that all the different available plot types may be listed by running plot-ms.py --list-plots.

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
plot1.png

plot1.png: our first plot: plot-ms.py 3C147_spw0.MS

plot2.0.png

plot2.0.png: second plot: plot-ms.py 3C147_spw0.MS DATA:XX DATA:YY

plot2.1.png

plot2.1.png: ...and page 2 of the above

plot3.png

plot3.png: plot-ms.py 3C147_spw0.MS DATA:XX.std

plot4.png

plot4.png: plot-ms.py 3C147_spw0.MS CORRECTED_DATA.mean:I

plot5.0.png

plot5.0.png: plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq

plot5.1.png

plot5.1.png: ...and page 2 of the above

plot6.0.png

plot6.0.png: plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq -F

plot6.1.png

plot6.1.png: ...and page 2 of the above

plot7.png

plot7.png: plot-ms.py 3C147_spw0.MS flags_XX

plot8.png

plot8.png: plot-ms.py 3C147_spw0.MS flags_XX --x-freq

diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/index.include.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/index.include.html new file mode 100644 index 0000000..7faa34b --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/index.include.html @@ -0,0 +1,82 @@ + +
+

Different kinds of plots

+

Logged on 02/27/2010 02:23:58 AM

+ + + +

Currently, the most versatile tool in the toolkit is the plot-ms.py script. To try it out, arm yourself with an MS. I will provide some examples based on the 3C147_spw0.MS, which you can download here:

+

http://www.astron.nl/meqwiki-data/users/oms/3C147-Calibration-Tutorial/3C147_spw0.MS.tgz

+

For starters, try

+

$ plot-ms.py -h

+

This just to let you know that there's many options, and I've tried to document them all. Don't bother reading them in detail, just type

+

$ plot-ms.py 3C147_spw0.MS

+

to make your first plot. This is a plot of the mean Stokes I (|XX|+|YY|) in the CORRECTED_DATA column, as a function of time, averaged over all frequencies. Note how data from individual baselines is stacked up and sorted by baseline length.

+

The plot is displayed in a window. Close the window to exit the script. To save the plot to a file, use the -o option:

+

$ plot-ms.py 3C147_spw0.MS -o plot1.png

+

Besides .png, you can also use .ps, .eps, .svg, and whatever else is supported by matplotlib -- see matplotlib documentation. For PostScript output, you can set paper size, plot size, etc -- see the help for details.

+

Now, this was a plot of mean Stokes I in CORRECTED_DATA. This is the default plot type made when nothing else is specified. To make different plots, you must specify them on the command line after the MS name, as "column:plottype". For example

+

$ plot-ms.py 3C147_spw0.MS DATA:XX DATA:YY

+

makes two plots, of mean |XX| and |YY| in the DATA column. Note the label at the top of the plot, which gives a bunch of details.

+

Instead of the mean |XX| in frequency, you can also plot standard deviations:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX.std

+

Other functions you can try are .min, .max and .sum. The default plot (plot 1) is equivalent to CORRECTED_DATA:I.mean. Finally, try:

+

$ plot-ms.py 3C147_spw0.MS CORRECTED_DATA.mean:I

+

and see if you can spot the difference. In plot 1, we plotted mean amplitudes, here we plot the amplitude of the (complex) mean visibility. This is an important distinction (for starters, the mean complex visibility always tends to be smaller, as it is a vector sum that tends to "wash out" due to the phase component.)! The positioning of ".mean" determines whether it applies to the column (i.e. visibilities), or to the amplitudes.

+

The previous plots were a function of time & averaged in frequency. To plot things as a function of frequency (averaged in time), use the --x-freq option. This will show you the time-frequency amplitude and phase bandpasses:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq

+

The gap in the plot is due to a flagged channel. To ignore flags (i.e. plot flagged data anyway), use the -F option. Note that this may impact the scale of your plot (which is chosen automatically based on the variation in the plotted data.)

+

Speaking of flags, another interesting plot type is the flag density plot. This shows the fraction of visibilities flagged per time slot or frequency channel. Use

+

$ plot-ms.py 3C147_spw0.MS flags_XX

+

to get a flag-density plot for XX (as a function of time.) Add the --x-freq option to plot it as a function of frequency. If you're flagging RFI, these kinds of plots are a nice way to see the bad timeslots and frequency channels.

+

Finally, note that all the different available plot types may be listed by running plot-ms.py --list-plots.

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
plot1.png

plot1.png: our first plot: plot-ms.py 3C147_spw0.MS

plot2.0.png

plot2.0.png: second plot: plot-ms.py 3C147_spw0.MS DATA:XX DATA:YY

plot2.1.png

plot2.1.png: ...and page 2 of the above

plot3.png

plot3.png: plot-ms.py 3C147_spw0.MS DATA:XX.std

plot4.png

plot4.png: plot-ms.py 3C147_spw0.MS CORRECTED_DATA.mean:I

plot5.0.png

plot5.0.png: plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq

plot5.1.png

plot5.1.png: ...and page 2 of the above

plot6.0.png

plot6.0.png: plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq -F

plot6.1.png

plot6.1.png: ...and page 2 of the above

plot7.png

plot7.png: plot-ms.py 3C147_spw0.MS flags_XX

plot8.png

plot8.png: plot-ms.py 3C147_spw0.MS flags_XX --x-freq

\ No newline at end of file diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot1.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot1.png new file mode 100644 index 0000000..7d0affc Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot1.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot1.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot1.png.purr-products/thumb.png new file mode 100644 index 0000000..3238661 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot1.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.0.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.0.png new file mode 100644 index 0000000..85b878e Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.0.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.0.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.0.png.purr-products/thumb.png new file mode 100644 index 0000000..09afee4 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.0.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.1.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.1.png new file mode 100644 index 0000000..2b44bb4 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.1.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.1.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.1.png.purr-products/thumb.png new file mode 100644 index 0000000..58a4c3f Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot2.1.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot3.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot3.png new file mode 100644 index 0000000..ad3df16 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot3.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot3.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot3.png.purr-products/thumb.png new file mode 100644 index 0000000..20e8476 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot3.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot4.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot4.png new file mode 100644 index 0000000..e67313c Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot4.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot4.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot4.png.purr-products/thumb.png new file mode 100644 index 0000000..5edf31f Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot4.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.0.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.0.png new file mode 100644 index 0000000..1220b99 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.0.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.0.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.0.png.purr-products/thumb.png new file mode 100644 index 0000000..9f206fa Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.0.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.1.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.1.png new file mode 100644 index 0000000..b594a61 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.1.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.1.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.1.png.purr-products/thumb.png new file mode 100644 index 0000000..3419d6f Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot5.1.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.0.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.0.png new file mode 100644 index 0000000..dd5d558 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.0.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.0.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.0.png.purr-products/thumb.png new file mode 100644 index 0000000..5cc6e1b Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.0.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.1.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.1.png new file mode 100644 index 0000000..8085778 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.1.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.1.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.1.png.purr-products/thumb.png new file mode 100644 index 0000000..17cd908 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot6.1.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot7.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot7.png new file mode 100644 index 0000000..8c8851c Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot7.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot7.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot7.png.purr-products/thumb.png new file mode 100644 index 0000000..b1e9664 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot7.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot8.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot8.png new file mode 100644 index 0000000..b659f49 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot8.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot8.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot8.png.purr-products/thumb.png new file mode 100644 index 0000000..40f30f7 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-022358/plot8.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/index.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/index.html new file mode 100644 index 0000000..7ac7028 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/index.html @@ -0,0 +1,60 @@ + + Data selection

Data selection

+

Logged on 02/27/2010 11:41:07 PM

+ + + +

Note the sawtooth pattern in the last flag density plot of the previous entry. This is due to this particular MS having every second channel flagged (in tribute to a long-gone imager bug.) We can use the -L option to select only the "good" channels.

+

$ plot-ms.py 3C147_spw0.MS flags_XX --x-freq -L 3~55:2

+

The -L option expects a single channel number, or a 'start~end' range ('end' is inclusive), or a 'start~end:step' specification like the above. If you're used to thinking in terms of Python slices (where N:M means from N up to and NOT including M), you can use the 'start:end' or 'start:end:step' syntax instead -- the first ':' indicates that 'end' is to be interpreted in the Pythonic "+1" sense. And, as you would expect, "10~" or "10:" means from 10 till the end, "10~:2" or "10::2" means from 10 till the end with a stepping of 2, and "~9" or ":10" means 0 to 9.

+

The -T option works the same way for timeslot numbers. To select the first 120 timeslots, try:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX -L 3:55:2 -T 0~119 # or 0:120, or just :120 for the Pythonic way

+

The next selection option is -I, to select interferometers. This uses the standard Meow.IfrSet object, and so has the same flexible syntax as the interferometer selector in all Siamese and Calico scripts. To get a summary of the syntax run:

+

$ plot-ms.py 3C147_spw0.MS -I help

+

In the examples below, we use -I "-<300* -C*" to select all baselines except those shorter than 300m and those with antenna C, "C* -CD" to select all baselines with C except C-D, and "F-M -7*" to select the standard WSRT set of fixed-movable baselines, but omit those with antenna 7.

+

Another data selection option is -D, to use with MSs that contain multiple DATA_DESC_IDs (these usually correspond to spectral windows.) The default is to use the first DDID in your MS (usually 0). You can also specify "-D all" or "-D 0,1,2" to plot data for multiple DDIDs (our test MS contains only DDID 0). Note that different DDIDs normally go on separate plot pages, but you can also ask the script to average them together. We'll learn about this later on in the tutorial.

+

Finally, the -Q option can be used to refine your data selection further through the use of TaQL queries. For example:

+

$ plot-ms.py 3C147_spw0.MS -Q "ANTENNA1==1"

+

is just another way of selecting all baselines where antenna 1 is the first element (equivalent to -I "1* -01"). TaQL is an SQL-like language for querying MSs. It is fully documented here: http://www.astron.nl/casacore/trunk/casacore/doc/notes/199.html.

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
plot9.png

plot9.png: plot-ms.py 3C147_spw0.MS flags_XX --x-freq -L 3~55:2

plot10.png

plot10.png: plot-ms.py 3C147_spw0.MS DATA:XX -L 3~55:2 -T 0~119

plot11.png

plot11.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''-<300 -C*''

plot12.png

plot12.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''C* -CD''

plot13.png

plot13.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''F-M -7*''

plot14.png

plot14.png: plot-ms.py 3C147_spw0.MS -Q ''ANTENNA1==1''

diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/index.include.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/index.include.html new file mode 100644 index 0000000..e90572c --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/index.include.html @@ -0,0 +1,49 @@ + +
+

Data selection

+

Logged on 02/27/2010 11:41:07 PM

+ + + +

Note the sawtooth pattern in the last flag density plot of the previous entry. This is due to this particular MS having every second channel flagged (in tribute to a long-gone imager bug.) We can use the -L option to select only the "good" channels.

+

$ plot-ms.py 3C147_spw0.MS flags_XX --x-freq -L 3~55:2

+

The -L option expects a single channel number, or a 'start~end' range ('end' is inclusive), or a 'start~end:step' specification like the above. If you're used to thinking in terms of Python slices (where N:M means from N up to and NOT including M), you can use the 'start:end' or 'start:end:step' syntax instead -- the first ':' indicates that 'end' is to be interpreted in the Pythonic "+1" sense. And, as you would expect, "10~" or "10:" means from 10 till the end, "10~:2" or "10::2" means from 10 till the end with a stepping of 2, and "~9" or ":10" means 0 to 9.

+

The -T option works the same way for timeslot numbers. To select the first 120 timeslots, try:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX -L 3:55:2 -T 0~119 # or 0:120, or just :120 for the Pythonic way

+

The next selection option is -I, to select interferometers. This uses the standard Meow.IfrSet object, and so has the same flexible syntax as the interferometer selector in all Siamese and Calico scripts. To get a summary of the syntax run:

+

$ plot-ms.py 3C147_spw0.MS -I help

+

In the examples below, we use -I "-<300* -C*" to select all baselines except those shorter than 300m and those with antenna C, "C* -CD" to select all baselines with C except C-D, and "F-M -7*" to select the standard WSRT set of fixed-movable baselines, but omit those with antenna 7.

+

Another data selection option is -D, to use with MSs that contain multiple DATA_DESC_IDs (these usually correspond to spectral windows.) The default is to use the first DDID in your MS (usually 0). You can also specify "-D all" or "-D 0,1,2" to plot data for multiple DDIDs (our test MS contains only DDID 0). Note that different DDIDs normally go on separate plot pages, but you can also ask the script to average them together. We'll learn about this later on in the tutorial.

+

Finally, the -Q option can be used to refine your data selection further through the use of TaQL queries. For example:

+

$ plot-ms.py 3C147_spw0.MS -Q "ANTENNA1==1"

+

is just another way of selecting all baselines where antenna 1 is the first element (equivalent to -I "1* -01"). TaQL is an SQL-like language for querying MSs. It is fully documented here: http://www.astron.nl/casacore/trunk/casacore/doc/notes/199.html.

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
plot9.png

plot9.png: plot-ms.py 3C147_spw0.MS flags_XX --x-freq -L 3~55:2

plot10.png

plot10.png: plot-ms.py 3C147_spw0.MS DATA:XX -L 3~55:2 -T 0~119

plot11.png

plot11.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''-<300 -C*''

plot12.png

plot12.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''C* -CD''

plot13.png

plot13.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''F-M -7*''

plot14.png

plot14.png: plot-ms.py 3C147_spw0.MS -Q ''ANTENNA1==1''

\ No newline at end of file diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot10.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot10.png new file mode 100644 index 0000000..989117f Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot10.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot10.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot10.png.purr-products/thumb.png new file mode 100644 index 0000000..8f1ae32 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot10.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot11.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot11.png new file mode 100644 index 0000000..0976d95 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot11.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot11.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot11.png.purr-products/thumb.png new file mode 100644 index 0000000..0811383 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot11.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot12.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot12.png new file mode 100644 index 0000000..56c8034 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot12.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot12.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot12.png.purr-products/thumb.png new file mode 100644 index 0000000..71fd100 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot12.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot13.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot13.png new file mode 100644 index 0000000..868a892 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot13.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot13.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot13.png.purr-products/thumb.png new file mode 100644 index 0000000..5a4140b Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot13.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot14.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot14.png new file mode 100644 index 0000000..31a9598 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot14.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot14.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot14.png.purr-products/thumb.png new file mode 100644 index 0000000..0adba12 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot14.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot9.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot9.png new file mode 100644 index 0000000..c43cca7 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot9.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot9.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot9.png.purr-products/thumb.png new file mode 100644 index 0000000..cd7bd31 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100227-234107/plot9.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/index.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/index.html new file mode 100644 index 0000000..a0c57e8 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/index.html @@ -0,0 +1,55 @@ + + Arranging plots

Arranging plots

+

Logged on 03/03/2010 12:46:53 AM

+ + + +

Normally, a command like

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*"

+

will produce two separate plot pages: one for |XX| and one for |YY|. The "-S" option can be used to stack the different plot types onto the same page. By default, only baselines are stacked, but the following command:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot

+

tells plot-ms to stack the plot types as well. To change the order in which the plots are stacked w.r.t.each other, you can give the -S option several times. For example, this will stack plots first by type, then by baseline:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot -S ifr

+

You can also put different baselines on different plot page by using the -P option. (Make sure you don't specify too many baselines in this case.) Try this:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "AB CD" -S plot -P ifr --size 20x10

+

(Note that we use the "--size" option to adjust the size of the plots.) Finally, the -A option can be used to average baselines together instead of stacking them:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -S plot -A ifr --size 20x10

+

And of course the same -P/-S/-A options can also apply to multiple DDIDs (i.e. spectral windows), if your MS happens to contain them.

+

This concludes our tutorial for the plot-ms tool. There are a few more obscure command-line options, you can use "plot-ms.py -h" to look them up.

+

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
plot15.png

plot15.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot

plot16.png

plot16.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot -S ifr

plot17.0.png

plot17.0.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "AB CD" -S plot -P ifr --size 20x10

plot17.1.png

plot17.1.png: page 2 of the above

plot18.png

plot18.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -S plot -A ifr --size 20x10

diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/index.include.html b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/index.include.html new file mode 100644 index 0000000..699c744 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/index.include.html @@ -0,0 +1,46 @@ + +
+

Arranging plots

+

Logged on 03/03/2010 12:46:53 AM

+ + + +

Normally, a command like

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*"

+

will produce two separate plot pages: one for |XX| and one for |YY|. The "-S" option can be used to stack the different plot types onto the same page. By default, only baselines are stacked, but the following command:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot

+

tells plot-ms to stack the plot types as well. To change the order in which the plots are stacked w.r.t.each other, you can give the -S option several times. For example, this will stack plots first by type, then by baseline:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot -S ifr

+

You can also put different baselines on different plot page by using the -P option. (Make sure you don't specify too many baselines in this case.) Try this:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "AB CD" -S plot -P ifr --size 20x10

+

(Note that we use the "--size" option to adjust the size of the plots.) Finally, the -A option can be used to average baselines together instead of stacking them:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -S plot -A ifr --size 20x10

+

And of course the same -P/-S/-A options can also apply to multiple DDIDs (i.e. spectral windows), if your MS happens to contain them.

+

This concludes our tutorial for the plot-ms tool. There are a few more obscure command-line options, you can use "plot-ms.py -h" to look them up.

+

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + +
plot15.png

plot15.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot

plot16.png

plot16.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot -S ifr

plot17.0.png

plot17.0.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "AB CD" -S plot -P ifr --size 20x10

plot17.1.png

plot17.1.png: page 2 of the above

plot18.png

plot18.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -S plot -A ifr --size 20x10

\ No newline at end of file diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot15.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot15.png new file mode 100644 index 0000000..a4aa93e Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot15.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot15.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot15.png.purr-products/thumb.png new file mode 100644 index 0000000..8f2bfc8 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot15.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot16.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot16.png new file mode 100644 index 0000000..28ee1ed Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot16.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot16.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot16.png.purr-products/thumb.png new file mode 100644 index 0000000..f4ad1f8 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot16.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.0.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.0.png new file mode 100644 index 0000000..f6a935b Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.0.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.0.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.0.png.purr-products/thumb.png new file mode 100644 index 0000000..2316c6a Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.0.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.1.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.1.png new file mode 100644 index 0000000..9a93cbd Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.1.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.1.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.1.png.purr-products/thumb.png new file mode 100644 index 0000000..f6e6a8b Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot17.1.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot18.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot18.png new file mode 100644 index 0000000..5859987 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot18.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot18.png.purr-products/thumb.png b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot18.png.purr-products/thumb.png new file mode 100644 index 0000000..eb61676 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/entry-20100303-004653/plot18.png.purr-products/thumb.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/ignore-20100227-030143/index.html b/tutorial/Owlcat-plotms-tutorial.purrlog/ignore-20100227-030143/index.html new file mode 100644 index 0000000..85a7a39 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/ignore-20100227-030143/index.html @@ -0,0 +1,21 @@ + + This is not a real log entry

This is not a real log entry

+

Logged on 02/27/2010 03:01:43 AM

+ + + +

+ + + + + + + + + + + + + + diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/ignorelist b/tutorial/Owlcat-plotms-tutorial.purrlog/ignorelist new file mode 100644 index 0000000..c2f1810 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/ignorelist @@ -0,0 +1,6 @@ +1267234749 ignore /home/oms/Waterhole/contrib/OMS/Owlcat/tutorial/plot3.1.png +1267234748 ignore /home/oms/Waterhole/contrib/OMS/Owlcat/tutorial/plot3.0.png +1267234607 ignore /home/oms/Waterhole/contrib/OMS/Owlcat/tutorial/plot2.1.png +1267234606 ignore /home/oms/Waterhole/contrib/OMS/Owlcat/tutorial/plot2.0.png +1267234847 ignore /home/oms/Waterhole/contrib/OMS/Owlcat/tutorial/plot3.png +1267235851 ignore /home/oms/Waterhole/contrib/OMS/Owlcat/tutorial/plot5.png diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/index.html b/tutorial/Owlcat-plotms-tutorial.purrlog/index.html new file mode 100644 index 0000000..81ccc43 --- /dev/null +++ b/tutorial/Owlcat-plotms-tutorial.purrlog/index.html @@ -0,0 +1,205 @@ + + + An Owlcat Tutorial + +

An Owlcat Tutorial

+ + +
+

Introducing Owlcat

+

Logged on 02/27/2010 02:12:41 AM

+ + + +

OWLCAT stands for Oleg's W--- Lightweight Calibration & Analysis Toolkit. The W can represent Wisely, Wantonly, Woefully, Weirdly, Wonderfully, Wickedly, Wobbly, Wankingly, Warp 9, or Whatever suits your mood. Owlcat is a bunch of scripts I wrote to help me reduce data. Of these, plot-ms.py, flag-ms.py and plot-parms.py are probably the most useful.

+

Owlcat is (currently) in my Waterhole, so if you want to use it, you must first check it out:

+

$ svn co svn://lofar9.astron.nl/var/svn/repos/trunk/Waterhole/contrib/OMS/Owlcat

+

This will create a directory called Owlcat. Now, make sure this directory is both in your PYTHONPATH and PATH.

+

Before you go any further, please make sure you've got Pylab/Matplotlib installed (Ubuntu package python-matplotlib), and have a copy of the MeqTrees Cattery scripts (svn co svn://lofar9.astron.nl/var/svn/repos/trunk/Frameworks/Cattery. No working MeqTrees required, just a copy of the Cattery python scripts.) Owlcat will not work without these.

+ + +

Data products

+ + + + + + +
owlcat-logo.jpg

owlcat-logo.jpg: The Owlcat logo

+
+

Different kinds of plots

+

Logged on 02/27/2010 02:23:58 AM

+ + + +

Currently, the most versatile tool in the toolkit is the plot-ms.py script. To try it out, arm yourself with an MS. I will provide some examples based on the 3C147_spw0.MS, which you can download here:

+

http://www.astron.nl/meqwiki-data/users/oms/3C147-Calibration-Tutorial/3C147_spw0.MS.tgz

+

For starters, try

+

$ plot-ms.py -h

+

This just to let you know that there's many options, and I've tried to document them all. Don't bother reading them in detail, just type

+

$ plot-ms.py 3C147_spw0.MS

+

to make your first plot. This is a plot of the mean Stokes I (|XX|+|YY|) in the CORRECTED_DATA column, as a function of time, averaged over all frequencies. Note how data from individual baselines is stacked up and sorted by baseline length.

+

The plot is displayed in a window. Close the window to exit the script. To save the plot to a file, use the -o option:

+

$ plot-ms.py 3C147_spw0.MS -o plot1.png

+

Besides .png, you can also use .ps, .eps, .svg, and whatever else is supported by matplotlib -- see matplotlib documentation. For PostScript output, you can set paper size, plot size, etc -- see the help for details.

+

Now, this was a plot of mean Stokes I in CORRECTED_DATA. This is the default plot type made when nothing else is specified. To make different plots, you must specify them on the command line after the MS name, as "column:plottype". For example

+

$ plot-ms.py 3C147_spw0.MS DATA:XX DATA:YY

+

makes two plots, of mean |XX| and |YY| in the DATA column. Note the label at the top of the plot, which gives a bunch of details.

+

Instead of the mean |XX| in frequency, you can also plot standard deviations:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX.std

+

Other functions you can try are .min, .max and .sum. The default plot (plot 1) is equivalent to CORRECTED_DATA:I.mean. Finally, try:

+

$ plot-ms.py 3C147_spw0.MS CORRECTED_DATA.mean:I

+

and see if you can spot the difference. In plot 1, we plotted mean amplitudes, here we plot the amplitude of the (complex) mean visibility. This is an important distinction (for starters, the mean complex visibility always tends to be smaller, as it is a vector sum that tends to "wash out" due to the phase component.)! The positioning of ".mean" determines whether it applies to the column (i.e. visibilities), or to the amplitudes.

+

The previous plots were a function of time & averaged in frequency. To plot things as a function of frequency (averaged in time), use the --x-freq option. This will show you the time-frequency amplitude and phase bandpasses:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq

+

The gap in the plot is due to a flagged channel. To ignore flags (i.e. plot flagged data anyway), use the -F option. Note that this may impact the scale of your plot (which is chosen automatically based on the variation in the plotted data.)

+

Speaking of flags, another interesting plot type is the flag density plot. This shows the fraction of visibilities flagged per time slot or frequency channel. Use

+

$ plot-ms.py 3C147_spw0.MS flags_XX

+

to get a flag-density plot for XX (as a function of time.) Add the --x-freq option to plot it as a function of frequency. If you're flagging RFI, these kinds of plots are a nice way to see the bad timeslots and frequency channels.

+

Finally, note that all the different available plot types may be listed by running plot-ms.py --list-plots.

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
plot1.png

plot1.png: our first plot: plot-ms.py 3C147_spw0.MS

plot2.0.png

plot2.0.png: second plot: plot-ms.py 3C147_spw0.MS DATA:XX DATA:YY

plot2.1.png

plot2.1.png: ...and page 2 of the above

plot3.png

plot3.png: plot-ms.py 3C147_spw0.MS DATA:XX.std

plot4.png

plot4.png: plot-ms.py 3C147_spw0.MS CORRECTED_DATA.mean:I

plot5.0.png

plot5.0.png: plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq

plot5.1.png

plot5.1.png: ...and page 2 of the above

plot6.0.png

plot6.0.png: plot-ms.py 3C147_spw0.MS DATA:XX XXphase --x-freq -F

plot6.1.png

plot6.1.png: ...and page 2 of the above

plot7.png

plot7.png: plot-ms.py 3C147_spw0.MS flags_XX

plot8.png

plot8.png: plot-ms.py 3C147_spw0.MS flags_XX --x-freq

+
+

Data selection

+

Logged on 02/27/2010 11:41:07 PM

+ + + +

Note the sawtooth pattern in the last flag density plot of the previous entry. This is due to this particular MS having every second channel flagged (in tribute to a long-gone imager bug.) We can use the -L option to select only the "good" channels.

+

$ plot-ms.py 3C147_spw0.MS flags_XX --x-freq -L 3~55:2

+

The -L option expects a single channel number, or a 'start~end' range ('end' is inclusive), or a 'start~end:step' specification like the above. If you're used to thinking in terms of Python slices (where N:M means from N up to and NOT including M), you can use the 'start:end' or 'start:end:step' syntax instead -- the first ':' indicates that 'end' is to be interpreted in the Pythonic "+1" sense. And, as you would expect, "10~" or "10:" means from 10 till the end, "10~:2" or "10::2" means from 10 till the end with a stepping of 2, and "~9" or ":10" means 0 to 9.

+

The -T option works the same way for timeslot numbers. To select the first 120 timeslots, try:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX -L 3:55:2 -T 0~119 # or 0:120, or just :120 for the Pythonic way

+

The next selection option is -I, to select interferometers. This uses the standard Meow.IfrSet object, and so has the same flexible syntax as the interferometer selector in all Siamese and Calico scripts. To get a summary of the syntax run:

+

$ plot-ms.py 3C147_spw0.MS -I help

+

In the examples below, we use -I "-<300* -C*" to select all baselines except those shorter than 300m and those with antenna C, "C* -CD" to select all baselines with C except C-D, and "F-M -7*" to select the standard WSRT set of fixed-movable baselines, but omit those with antenna 7.

+

Another data selection option is -D, to use with MSs that contain multiple DATA_DESC_IDs (these usually correspond to spectral windows.) The default is to use the first DDID in your MS (usually 0). You can also specify "-D all" or "-D 0,1,2" to plot data for multiple DDIDs (our test MS contains only DDID 0). Note that different DDIDs normally go on separate plot pages, but you can also ask the script to average them together. We'll learn about this later on in the tutorial.

+

Finally, the -Q option can be used to refine your data selection further through the use of TaQL queries. For example:

+

$ plot-ms.py 3C147_spw0.MS -Q "ANTENNA1==1"

+

is just another way of selecting all baselines where antenna 1 is the first element (equivalent to -I "1* -01"). TaQL is an SQL-like language for querying MSs. It is fully documented here: http://www.astron.nl/casacore/trunk/casacore/doc/notes/199.html.

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
plot9.png

plot9.png: plot-ms.py 3C147_spw0.MS flags_XX --x-freq -L 3~55:2

plot10.png

plot10.png: plot-ms.py 3C147_spw0.MS DATA:XX -L 3~55:2 -T 0~119

plot11.png

plot11.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''-<300 -C*''

plot12.png

plot12.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''C* -CD''

plot13.png

plot13.png: plot-ms.py 3C147_spw0.MS DATA:XX -I ''F-M -7*''

plot14.png

plot14.png: plot-ms.py 3C147_spw0.MS -Q ''ANTENNA1==1''

+
+

Arranging plots

+

Logged on 03/03/2010 12:46:53 AM

+ + + +

Normally, a command like

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*"

+

will produce two separate plot pages: one for |XX| and one for |YY|. The "-S" option can be used to stack the different plot types onto the same page. By default, only baselines are stacked, but the following command:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot

+

tells plot-ms to stack the plot types as well. To change the order in which the plots are stacked w.r.t.each other, you can give the -S option several times. For example, this will stack plots first by type, then by baseline:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot -S ifr

+

You can also put different baselines on different plot page by using the -P option. (Make sure you don't specify too many baselines in this case.) Try this:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -I "AB CD" -S plot -P ifr --size 20x10

+

(Note that we use the "--size" option to adjust the size of the plots.) Finally, the -A option can be used to average baselines together instead of stacking them:

+

$ plot-ms.py 3C147_spw0.MS DATA:XX YY -S plot -A ifr --size 20x10

+

And of course the same -P/-S/-A options can also apply to multiple DDIDs (i.e. spectral windows), if your MS happens to contain them.

+

This concludes our tutorial for the plot-ms tool. There are a few more obscure command-line options, you can use "plot-ms.py -h" to look them up.

+

+ + +

Data products

+ + + + + + + + + + + + + + + + + + + + + + +
plot15.png

plot15.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot

plot16.png

plot16.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "C*" -S plot -S ifr

plot17.0.png

plot17.0.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -I "AB CD" -S plot -P ifr --size 20x10

plot17.1.png

plot17.1.png: page 2 of the above

plot18.png

plot18.png: plot-ms.py 3C147_spw0.MS DATA:XX YY -S plot -A ifr --size 20x10

+
This log was generated + by PURR version 1.1.
+ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/purr16.png b/tutorial/Owlcat-plotms-tutorial.purrlog/purr16.png new file mode 100644 index 0000000..692dd0d Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/purr16.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/purr24.png b/tutorial/Owlcat-plotms-tutorial.purrlog/purr24.png new file mode 100644 index 0000000..40d2110 Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/purr24.png differ diff --git a/tutorial/Owlcat-plotms-tutorial.purrlog/purr32.png b/tutorial/Owlcat-plotms-tutorial.purrlog/purr32.png new file mode 100644 index 0000000..af81e0b Binary files /dev/null and b/tutorial/Owlcat-plotms-tutorial.purrlog/purr32.png differ