source: branches/Release12/python/scantable.py @ 2569

Last change on this file since 2569 was 800, checked in by mar637, 18 years ago

Request: added channel based flagging
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 63.9 KB
Line 
1from asap._asap import sdtable
2from asap import rcParams
3from asap import print_log
4from numarray import ones,zeros
5import sys
6
7class scantable(sdtable):
8    """
9        The ASAP container for scans
10    """
11
12    def __init__(self, filename, average=None, unit=None):
13        """
14        Create a scantable from a saved one or make a reference
15        Parameters:
16            filename:    the name of an asap table on disk
17                         or
18                         the name of a rpfits/sdfits/ms file
19                         (integrations within scans are auto averaged
20                         and the whole file is read)
21                         or
22                         [advanced] a reference to an existing
23                         scantable
24            average:     average all integrations withinb a scan on read.
25                         The default (True) is taken from .asaprc.
26            unit:         brightness unit; must be consistent with K or Jy.
27                         Over-rides the default selected by the reader
28                         (input rpfits/sdfits/ms) or replaces the value
29                         in existing scantables
30        """
31        if average is None or type(average) is not bool:
32            average = rcParams['scantable.autoaverage']
33
34        varlist = vars()
35        self._p = None
36        from asap import asaplog
37        if isinstance(filename,sdtable):
38            sdtable.__init__(self, filename)
39            if unit is not None:
40                self.set_fluxunit(unit)
41        else:
42            import os.path
43            if not os.path.exists(filename):
44                s = "File '%s' not found." % (filename)
45                if rcParams['verbose']:
46                    asaplog.push(s)
47                    print asaplog.pop().strip()
48                    return
49                raise IOError(s)
50            filename = os.path.expandvars(filename)
51            if os.path.isdir(filename):
52                # crude check if asap table
53                if os.path.exists(filename+'/table.info'):
54                    sdtable.__init__(self, filename)
55                    if unit is not None:
56                        self.set_fluxunit(unit)
57                else:
58                    msg = "The given file '%s'is not a valid asap table." % (filename)
59                    if rcParams['verbose']:
60                        print msg
61                        return
62                    else:
63                        raise IOError(msg)
64            else:
65                from asap._asap import sdreader
66                ifSel = -1
67                beamSel = -1
68                r = sdreader()
69                r._open(filename,ifSel,beamSel)
70                asaplog.push('Importing data...')
71                print_log()
72                r._read([-1])
73                tbl = r._getdata()
74                if unit is not None:
75                    tbl.set_fluxunit(unit)
76                if average:
77                    from asap._asap import average as _av
78                    asaplog.push('Auto averaging integrations...')
79                    print_log()
80                    tbl2 = _av((tbl,),(),True,'none')
81                    sdtable.__init__(self,tbl2)
82                    del tbl2
83                else:
84                    sdtable.__init__(self,tbl)
85                del r,tbl
86                self._add_history("scantable", varlist)
87        print_log()
88
89    def save(self, name=None, format=None, stokes=False, overwrite=False):
90        """
91        Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
92        Image FITS or MS2 format.
93        Parameters:
94            name:        the name of the outputfile. For format="FITS" this
95                         is the directory file name into which all the files
96                         will be written (default is 'asap_FITS'). For format
97                         "ASCII" this is the root file name (data in 'name'.txt
98                         and header in 'name'_header.txt)
99            format:      an optional file format. Default is ASAP.
100                         Allowed are - 'ASAP' (save as ASAP [aips++] Table),
101                                       'SDFITS' (save as SDFITS file)
102                                       'FITS' (saves each row as a FITS Image)
103                                       'ASCII' (saves as ascii text file)
104                                       'MS2' (saves as an aips++
105                                              MeasurementSet V2)
106            stokes:      Convert to Stokes parameters (only available
107                         currently with FITS and ASCII formats.
108                         Default is False.
109            overwrite:   If the file should be overwritten if it exists.
110                         The default False is to return with warning
111                         without writing the output. USE WITH CARE.
112        Example:
113            scan.save('myscan.asap')
114            scan.save('myscan.sdfits','SDFITS')
115        """
116        from os import path
117        if format is None: format = rcParams['scantable.save']
118        suffix = '.'+format.lower()
119        if name is None or name =="":
120            name = 'scantable'+suffix
121            from asap import asaplog
122            msg = "No filename given. Using default name %s..." % name
123            asaplog.push(msg)
124        name = path.expandvars(name)
125        if path.isfile(name) or path.isdir(name):
126            if not overwrite:
127                msg = "File %s exists." % name
128                if rcParams['verbose']:
129                    print msg
130                    return
131                else:
132                    raise IOError(msg)
133        format2 = format.upper()
134        if format2 == 'ASAP':
135            self._save(name)
136        else:
137            from asap._asap import sdwriter as _sw
138            w = _sw(format2)
139            w.write(self, name, stokes)
140
141        print_log()
142        return
143
144    def copy(self):
145        """
146        Return a copy of this scantable.
147        Parameters:
148            none
149        Example:
150            copiedscan = scan.copy()
151        """
152        sd = scantable(sdtable._copy(self))
153        return sd
154
155    def get_scan(self, scanid=None):
156        """
157        Return a specific scan (by scanno) or collection of scans (by
158        source name) in a new scantable.
159        Parameters:
160            scanid:    a (list of) scanno or a source name, unix-style
161                       patterns are accepted for source name matching, e.g.
162                       '*_R' gets all 'ref scans
163        Example:
164            # get all scans containing the source '323p459'
165            newscan = scan.get_scan('323p459')
166            # get all 'off' scans
167            refscans = scan.get_scan('*_R')
168            # get a susbset of scans by scanno (as listed in scan.summary())
169            newscan = scan.get_scan([0,2,7,10])
170        """
171        if scanid is None:
172            if rcParams['verbose']:
173                print "Please specify a scan no or name to retrieve from the scantable"
174                return
175            else:
176                raise RuntimeError("No scan given")
177
178        try:
179            if type(scanid) is str:
180                s = sdtable._getsource(self,scanid)
181                return scantable(s)
182            elif type(scanid) is int:
183                s = sdtable._getscan(self,[scanid])
184                return scantable(s)
185            elif type(scanid) is list:
186                s = sdtable._getscan(self,scanid)
187                return scantable(s)
188            else:
189                msg = "Illegal scanid type, use 'int' or 'list' if ints."
190                if rcParams['verbose']:
191                    print msg
192                else:
193                    raise TypeError(msg)
194        except RuntimeError:
195            if rcParams['verbose']: print "Couldn't find any match."
196            else: raise
197
198    def __str__(self):
199        return sdtable._summary(self,True)
200
201    def summary(self,filename=None, verbose=None):
202        """
203        Print a summary of the contents of this scantable.
204        Parameters:
205            filename:    the name of a file to write the putput to
206                         Default - no file output
207            verbose:     print extra info such as the frequency table
208                         The default (False) is taken from .asaprc
209        """
210        info = sdtable._summary(self, verbose)
211        if verbose is None: verbose = rcParams['scantable.verbosesummary']
212        if filename is not None:
213            if filename is "":
214                filename = 'scantable_summary.txt'
215            from os.path import expandvars, isdir
216            filename = expandvars(filename)
217            if not isdir(filename):
218                data = open(filename, 'w')
219                data.write(info)
220                data.close()
221            else:
222                msg = "Illegal file name '%s'." % (filename)
223                if rcParams['verbose']:
224                    print msg
225                else:
226                    raise IOError(msg)
227        if rcParams['verbose']:
228            try:
229                from IPython.genutils import page as pager
230            except ImportError:
231                from pydoc import pager
232            pager(info)
233        else:
234            return info
235
236    def set_cursor(self, beam=0, IF=0, pol=0):
237        """
238        Set the spectrum for individual operations.
239        Parameters:
240            beam, IF, pol:    a number
241        Example:
242            scan.set_cursor(0,0,1)
243            pol1sig = scan.stats(all=False) # returns std dev for beam=0
244                                            # if=0, pol=1
245        """
246        varlist = vars()
247        self.setbeam(beam)
248        self.setpol(pol)
249        self.setif(IF)
250        self._add_history("set_cursor",varlist)
251        return
252
253    def get_cursor(self):
254        """
255        Return/print a the current 'cursor' into the Beam/IF/Pol cube.
256        Parameters:
257            none
258        Returns:
259            a list of values (currentBeam,currentIF,currentPol)
260        Example:
261            none
262        """
263        i = self.getbeam()
264        j = self.getif()
265        k = self.getpol()
266        from asap import asaplog
267        out = "--------------------------------------------------\n"
268        out += " Cursor position\n"
269        out += "--------------------------------------------------\n"
270        out += 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
271        asaplog.push(out)
272        print_log()
273        return i,j,k
274
275    def stats(self, stat='stddev', mask=None, allaxes=None):
276        """
277        Determine the specified statistic of the current beam/if/pol
278        Takes a 'mask' as an optional parameter to specify which
279        channels should be excluded.
280        Parameters:
281            stat:    'min', 'max', 'sumsq', 'sum', 'mean'
282                     'var', 'stddev', 'avdev', 'rms', 'median'
283            mask:    an optional mask specifying where the statistic
284                     should be determined.
285            allaxes: if True apply to all spectra. Otherwise
286                     apply only to the selected (beam/pol/if)spectra only.
287                     The default is taken from .asaprc (True if none)
288        Example:
289            scan.set_unit('channel')
290            msk = scan.create_mask([100,200],[500,600])
291            scan.stats(stat='mean', mask=m)
292        """
293        if allaxes is None: allaxes = rcParams['scantable.allaxes']
294        from asap._asap import stats as _stats
295        from numarray import array,zeros,Float
296        if mask == None:
297            mask = ones(self.nchan())
298        axes = ['Beam','IF','Pol','Time']
299
300        beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
301        if allaxes:
302            n = self.nbeam()*self.nif()*self.npol()*self.nrow()
303            shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
304            arr = array(zeros(n),shape=shp,type=Float)
305
306            for i in range(self.nbeam()):
307                self.setbeam(i)
308                for j in range(self.nif()):
309                    self.setif(j)
310                    for k in range(self.npol()):
311                        self.setpol(k)
312                        arr[i,j,k,:] = _stats(self,mask,stat,-1)
313            retval = {'axes': axes, 'data': arr, 'cursor':None}
314            tm = [self._gettime(val) for val in range(self.nrow())]
315            if rcParams['verbose']:
316                self._print_values(retval,stat,tm)
317            self.setbeam(beamSel)
318            self.setif(IFSel)
319            self.setpol(polSel)
320            return retval
321
322        else:
323            statval = _stats(self,mask,stat,-1)
324            out = ''
325            for l in range(self.nrow()):
326                tm = self._gettime(l)
327                out += 'Time[%s]:\n' % (tm)
328                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (beamSel)
329                if self.nif() > 1: out +=  ' IF[%d] ' % (IFSel)
330                if self.npol() > 1: out +=  ' Pol[%d] ' % (polSel)
331                out += '= %3.3f\n' % (statval[l])
332                out +=  "--------------------------------------------------\n"
333
334            if rcParams['verbose']:
335                print "--------------------------------------------------"
336                print " ",stat
337                print "--------------------------------------------------"
338                print out
339            retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
340            return retval
341
342    def stddev(self,mask=None, allaxes=None):
343        """
344        Determine the standard deviation of the current beam/if/pol
345        Takes a 'mask' as an optional parameter to specify which
346        channels should be excluded.
347        Parameters:
348            mask:    an optional mask specifying where the standard
349                     deviation should be determined.
350            allaxes: optional flag to show all or a cursor selected
351                     spectrum of Beam/IF/Pol. Default is all or taken
352                     from .asaprc
353
354        Example:
355            scan.set_unit('channel')
356            msk = scan.create_mask([100,200],[500,600])
357            scan.stddev(mask=m)
358        """
359        if allaxes is None: allaxes = rcParams['scantable.allaxes']
360        return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
361
362    def get_tsys(self, allaxes=None):
363        """
364        Return the System temperatures.
365        Parameters:
366           allaxes:     if True apply to all spectra. Otherwise
367                        apply only to the selected (beam/pol/if)spectra only.
368                        The default is taken from .asaprc (True if none)
369        Returns:
370            a list of Tsys values.
371        """
372        if allaxes is None: allaxes = rcParams['scantable.allaxes']
373        from numarray import array,zeros,Float
374        axes = ['Beam','IF','Pol','Time']
375
376        if allaxes:
377            n = self.nbeam()*self.nif()*self.npol()*self.nrow()
378            shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
379            arr = array(zeros(n),shape=shp,type=Float)
380
381            for i in range(self.nbeam()):
382                self.setbeam(i)
383                for j in range(self.nif()):
384                    self.setif(j)
385                    for k in range(self.npol()):
386                        self.setpol(k)
387                        arr[i,j,k,:] = self._gettsys()
388            retval = {'axes': axes, 'data': arr, 'cursor':None}
389            tm = [self._gettime(val) for val in range(self.nrow())]
390            if rcParams['verbose']:
391                self._print_values(retval,'Tsys',tm)
392            return retval
393
394        else:
395            i,j,k = (self.getbeam(),self.getif(),self.getpol())
396            statval = self._gettsys()
397            out = ''
398            for l in range(self.nrow()):
399                tm = self._gettime(l)
400                out += 'Time[%s]:\n' % (tm)
401                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
402                if self.nif() > 1: out +=  ' IF[%d] ' % (j)
403                if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
404                out += '= %3.3f\n' % (statval[l])
405                out +=  "--------------------------------------------------\n"
406
407            if rcParams['verbose']:
408                print "--------------------------------------------------"
409                print " TSys"
410                print "--------------------------------------------------"
411                print out
412            retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
413            return retval
414
415    def get_time(self, row=-1):
416        """
417        Get a list of time stamps for the observations.
418        Return a string for each integration in the scantable.
419        Parameters:
420            row:    row no of integration. Default -1 return all rows
421        Example:
422            none
423        """
424        out = []
425        if row == -1:
426            for i in range(self.nrow()):
427                out.append(self._gettime(i))
428            return out
429        else:
430            if row < self.nrow():
431                return self._gettime(row)
432
433    def get_sourcename(self, row=-1):
434        """
435        Get a list source names for the observations.
436        Return a string for each integration in the scantable.
437        Parameters:
438            row:    row no of integration. Default -1 return all rows
439        Example:
440            none
441        """
442        out = []
443        if row == -1:
444            return [self._getsourcename(i) for i in range(self.nrow())]
445        else:
446            if  0 <= row < self.nrow():
447                return self._getsourcename(row)
448
449    def get_elevation(self, row=-1):
450        """
451        Get a list of elevations for the observations.
452        Return a float for each integration in the scantable.
453        Parameters:
454            row:    row no of integration. Default -1 return all rows
455        Example:
456            none
457        """
458        out = []
459        if row == -1:
460            return [self._getelevation(i) for i in range(self.nrow())]
461        else:
462            if  0 <= row < self.nrow():
463                return self._getelevation(row)
464
465    def get_azimuth(self, row=-1):
466        """
467        Get a list of azimuths for the observations.
468        Return a float for each integration in the scantable.
469        Parameters:
470            row:    row no of integration. Default -1 return all rows
471        Example:
472            none
473        """
474        out = []
475        if row == -1:
476            return [self._getazimuth(i) for i in range(self.nrow())]
477        else:
478            if  0 <= row < self.nrow():
479                return self._getazimuth(row)
480
481    def get_parangle(self, row=-1):
482        """
483        Get a list of parallactic angles for the observations.
484        Return a float for each integration in the scantable.
485        Parameters:
486            row:    row no of integration. Default -1 return all rows
487        Example:
488            none
489        """
490        out = []
491        if row == -1:
492            return [self._getparangle(i) for i in range(self.nrow())]
493        else:
494            if  0 <= row < self.nrow():
495                return self._getparangle(row)
496
497    def set_unit(self, unit='channel'):
498        """
499        Set the unit for all following operations on this scantable
500        Parameters:
501            unit:    optional unit, default is 'channel'
502                     one of '*Hz','km/s','channel', ''
503        """
504        varlist = vars()
505        if unit in ['','pixel', 'channel']:
506            unit = ''
507        inf = list(self._getcoordinfo())
508        inf[0] = unit
509        self._setcoordinfo(inf)
510        if self._p: self.plot()
511        self._add_history("set_unit",varlist)
512
513    def set_instrument(self, instr):
514        """
515        Set the instrument for subsequent processing
516        Parameters:
517            instr:    Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
518                      'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
519        """
520        self._setInstrument(instr)
521        self._add_history("set_instument",vars())
522        print_log()
523
524    def set_doppler(self, doppler='RADIO'):
525        """
526        Set the doppler for all following operations on this scantable.
527        Parameters:
528            doppler:    One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
529        """
530        varlist = vars()
531        inf = list(self._getcoordinfo())
532        inf[2] = doppler
533        self._setcoordinfo(inf)
534        if self._p: self.plot()
535        self._add_history("set_doppler",vars())
536        print_log()
537
538    def set_freqframe(self, frame=None):
539        """
540        Set the frame type of the Spectral Axis.
541        Parameters:
542            frame:   an optional frame type, default 'LSRK'. Valid frames are:
543                     'REST','TOPO','LSRD','LSRK','BARY',
544                     'GEO','GALACTO','LGROUP','CMB'
545        Examples:
546            scan.set_freqframe('BARY')
547        """
548        if frame is None: frame = rcParams['scantable.freqframe']
549        varlist = vars()
550        valid = ['REST','TOPO','LSRD','LSRK','BARY', \
551                   'GEO','GALACTO','LGROUP','CMB']
552
553        if frame in valid:
554            inf = list(self._getcoordinfo())
555            inf[1] = frame
556            self._setcoordinfo(inf)
557            self._add_history("set_freqframe",varlist)
558        else:
559            msg  = "Please specify a valid freq type. Valid types are:\n",valid
560            if rcParams['verbose']:
561                print msg
562            else:
563                raise TypeError(msg)
564        print_log()
565
566    def get_unit(self):
567        """
568        Get the default unit set in this scantable
569        Parameters:
570        Returns:
571            A unit string
572        """
573        inf = self._getcoordinfo()
574        unit = inf[0]
575        if unit == '': unit = 'channel'
576        return unit
577
578    def get_abcissa(self, rowno=0):
579        """
580        Get the abcissa in the current coordinate setup for the currently
581        selected Beam/IF/Pol
582        Parameters:
583            rowno:    an optional row number in the scantable. Default is the
584                      first row, i.e. rowno=0
585        Returns:
586            The abcissa values and it's format string (as a dictionary)
587        """
588        abc = self._getabcissa(rowno)
589        lbl = self._getabcissalabel(rowno)
590        print_log()
591        return abc, lbl
592
593    def create_mask(self, *args, **kwargs):
594        """
595        Compute and return a mask based on [min,max] windows.
596        The specified windows are to be INCLUDED, when the mask is
597        applied.
598        Parameters:
599            [min,max],[min2,max2],...
600                Pairs of start/end points specifying the regions
601                to be masked
602            invert:     optional argument. If specified as True,
603                        return an inverted mask, i.e. the regions
604                        specified are EXCLUDED
605            row:        create the mask using the specified row for
606                        unit conversions, default is row=0
607                        only necessary if frequency varies over rows.
608        Example:
609            scan.set_unit('channel')
610
611            a)
612            msk = scan.set_mask([400,500],[800,900])
613            # masks everything outside 400 and 500
614            # and 800 and 900 in the unit 'channel'
615
616            b)
617            msk = scan.set_mask([400,500],[800,900], invert=True)
618            # masks the regions between 400 and 500
619            # and 800 and 900 in the unit 'channel'
620
621        """
622        row = 0
623        if kwargs.has_key("row"):
624            row = kwargs.get("row")
625        data = self._getabcissa(row)
626        u = self._getcoordinfo()[0]
627        if rcParams['verbose']:
628            if u == "": u = "channel"
629            from asap import asaplog
630            msg = "The current mask window unit is %s" % u
631            asaplog.push(msg)
632        n = self.nchan()
633        msk = zeros(n)
634        # test if args is a 'list' or a 'normal *args - UGLY!!!
635
636        ws = (isinstance(args[-1][-1],int) or isinstance(args[-1][-1],float)) and args or args[0]
637        for window in ws:
638            if (len(window) != 2 or window[0] > window[1] ):
639                raise TypeError("A window needs to be defined as [min,max]")
640            for i in range(n):
641                if data[i] >= window[0] and data[i] < window[1]:
642                    msk[i] = 1
643        if kwargs.has_key('invert'):
644            if kwargs.get('invert'):
645                from numarray import logical_not
646                msk = logical_not(msk)
647        print_log()
648        return msk
649
650    def get_restfreqs(self):
651        """
652        Get the restfrequency(s) stored in this scantable.
653        The return value(s) are always of unit 'Hz'
654        Parameters:
655            none
656        Returns:
657            a list of doubles
658        """
659        return list(self._getrestfreqs())
660
661    def lines(self):
662        """
663        Print the list of known spectral lines
664        """
665        l = sdtable._lines(self)
666        if rcParams['verbose']:
667            print l
668        else:
669            return l
670
671    def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
672                      theif=None):
673        """
674        Select the restfrequency for the specified source and IF OR
675        replace for all IFs.  If the 'freqs' argument holds a scalar,
676        then that rest frequency will be applied to the selected
677        data (and added to the list of available rest frequencies).
678        In this way, you can set a rest frequency for each
679        source and IF combination.   If the 'freqs' argument holds
680        a vector, then it MUST be of length the number of IFs
681        (and the available restfrequencies will be replaced by
682        this vector).  In this case, *all* data ('source' and
683        'theif' are ignored) have the restfrequency set per IF according
684        to the corresponding value you give in the 'freqs' vector.
685        E.g. 'freqs=[1e9,2e9]'  would mean IF 0 gets restfreq 1e9 and
686        IF 1 gets restfreq 2e9.
687
688        You can also specify the frequencies via known line names
689        in the argument 'lines'.  Use 'freqs' or 'lines'.  'freqs'
690        takes precedence. See the list of known names via function
691        scantable.lines()
692        Parameters:
693            freqs:   list of rest frequencies
694            unit:    unit for rest frequency (default 'Hz')
695            lines:   list of known spectral lines names (alternative to freqs).
696                     See possible list via scantable.lines()
697            source:  Source name (blank means all)
698            theif:   IF (-1 means all)
699        Example:
700            scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
701            scan.set_restfreqs(freqs=[1.4e9,1.67e9])
702        """
703        varlist = vars()
704        if source is None:
705            source = ""
706        if theif is None:
707            theif = -1
708        t = type(freqs)
709        if t is int or t is float:
710           freqs = [freqs]
711        if freqs is None:
712           freqs = []
713        t = type(lines)
714        if t is str:
715           lines = [lines]
716        if lines is None:
717           lines = []
718        sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
719        self._add_history("set_restfreqs", varlist)
720
721
722
723    def flag_spectrum(self, flags, row=-1, allaxes=None):
724        """
725        This flags a selected spectrum in the scan 'for good'.
726        USE WITH CARE - not reversible.
727        Use masks for non-permanent exclusion of channels.
728        Parameters:
729            flags:   a mask created with scantable.create_mask
730            row:     the number of the row (intergration) to flag,
731                     the deafult is -1, which is ALL rows.
732            allaxes: if True apply to all spectra. Otherwise
733                     apply only to the selected (beam/pol/if)spectra only
734                     The default is taken from .asaprc (True if none)
735        Example:
736            scan.flag_spectrum() # flags everything, not very useful
737            msk = scan.create_mask([511,513])
738            scan,set_cursor(IF=1)
739            scan.flag_spectrum(msk, allaxes=False)
740            # flags the central three channels for all rows in the table,
741            # using the selected IF only (e.g. a birdie)
742        """
743        if allaxes is None: allaxes = rcParams['scantable.allaxes']
744        if len(flags) != self.nchan():
745            msg  = "Length of mask invalid.",valid
746            if rcParams['verbose']:
747                print msg
748                return
749            else:
750                raise AssertionError(msg)
751        rows = range(self.nrow())
752        if row > -1: rows = [row]
753        for r in rows:
754            if allaxes:
755                beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
756                for i in range(self.nbeam()):
757                    self.setbeam(i)
758                    for j in range(self.nif()):
759                        self.setif(j)
760                        for k in range(self.npol()):
761                            self.setpol(k)
762                            sdtable._flagspectrum(self, flags, r)
763                self.setbeam(beamSel)
764                self.setif(IFSel)
765                self.setpol(polSel)
766            else:
767                sdtable._flagspectrum(self, flags, r)
768        self._add_history("flag_spectrum", vars())
769        return
770
771    def _print_values(self, dat, label='', timestamps=[]):
772        d = dat['data']
773        a = dat['axes']
774        shp = d.getshape()
775        out = ''
776        for i in range(shp[3]):
777            out += '%s [%s]:\n' % (a[3],timestamps[i])
778            t = d[:,:,:,i]
779            for j in range(shp[0]):
780                if shp[0] > 1: out +=  ' %s[%d] ' % (a[0],j)
781                for k in range(shp[1]):
782                    if shp[1] > 1: out +=  ' %s[%d] ' % (a[1],k)
783                    for l in range(shp[2]):
784                        if shp[2] > 1: out +=  ' %s[%d] ' % (a[2],l)
785                        out += '= %3.3f\n' % (t[j,k,l])
786            out += "-"*80
787            out += "\n"
788        print "-"*80
789        print " ", label
790        print "-"*80
791        print out
792
793    def history(self):
794        hist = list(self._gethistory())
795        out = "-"*80
796        for h in hist:
797            if h.startswith("---"):
798                out += "\n"+h
799            else:
800                items = h.split("##")
801                date = items[0]
802                func = items[1]
803                items = items[2:]
804                out += "\n"+date+"\n"
805                out += "Function: %s\n  Parameters:" % (func)
806                for i in items:
807                    s = i.split("=")
808                    out += "\n   %s = %s" % (s[0],s[1])
809                out += "\n"+"-"*80
810        try:
811            from IPython.genutils import page as pager
812        except ImportError:
813            from pydoc import pager
814        pager(out)
815        return
816
817    #
818    # Maths business
819    #
820
821    def average_time(self, mask=None, scanav=False, weight='tint'):
822        """
823        Return the (time) average of a scan, or apply it 'insitu'.
824        Note:
825            in channels only
826            The cursor of the output scan is set to 0.
827        Parameters:
828            one scan or comma separated  scans
829            mask:     an optional mask (only used for 'var' and 'tsys'
830                      weighting)
831            scanav:   True averages each scan separately
832                      False (default) averages all scans together,
833            weight:   Weighting scheme. 'none', 'var' (1/var(spec)
834                      weighted), 'tsys' (1/Tsys**2 weighted), 'tint'
835                      (integration time weighted) or 'tintsys' (Tint/Tsys**2).
836                      The default is 'tint'
837        Example:
838            # time average the scantable without using a mask
839            newscan = scan.average_time()
840        """
841        varlist = vars()
842        if weight is None: weight = 'tint'
843        if mask is None: mask = ()
844        from asap._asap import average as _av
845        s = scantable(_av((self,), mask, scanav, weight))
846        s._add_history("average_time",varlist)
847        print_log()
848        return s
849
850    def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
851                     allaxes=None):
852        """
853        Return a scan where all spectra are converted to either
854        Jansky or Kelvin depending upon the flux units of the scan table.
855        By default the function tries to look the values up internally.
856        If it can't find them (or if you want to over-ride), you must
857        specify EITHER jyperk OR eta (and D which it will try to look up
858        also if you don't set it). jyperk takes precedence if you set both.
859        Parameters:
860            jyperk:      the Jy / K conversion factor
861            eta:         the aperture efficiency
862            d:           the geomtric diameter (metres)
863            insitu:      if False a new scantable is returned.
864                         Otherwise, the scaling is done in-situ
865                         The default is taken from .asaprc (False)
866            allaxes:     if True apply to all spectra. Otherwise
867                         apply only to the selected (beam/pol/if)spectra only
868                         The default is taken from .asaprc (True if none)
869        """
870        if allaxes is None: allaxes = rcParams['scantable.allaxes']
871        if insitu is None: insitu = rcParams['insitu']
872        varlist = vars()
873        if jyperk is None: jyperk = -1.0
874        if d is None: d = -1.0
875        if eta is None: eta = -1.0
876        if not insitu:
877            from asap._asap import convertflux as _convert
878            s = scantable(_convert(self, d, eta, jyperk, allaxes))
879            s._add_history("convert_flux", varlist)
880            print_log()
881            return s
882        else:
883            from asap._asap import convertflux_insitu as _convert
884            _convert(self, d, eta, jyperk, allaxes)
885            self._add_history("convert_flux", varlist)
886            print_log()
887            return
888
889    def gain_el(self, poly=None, filename="", method="linear",
890                insitu=None, allaxes=None):
891        """
892        Return a scan after applying a gain-elevation correction.
893        The correction can be made via either a polynomial or a
894        table-based interpolation (and extrapolation if necessary).
895        You specify polynomial coefficients, an ascii table or neither.
896        If you specify neither, then a polynomial correction will be made
897        with built in coefficients known for certain telescopes (an error
898        will occur if the instrument is not known).
899        The data and Tsys are *divided* by the scaling factors.
900        Parameters:
901            poly:        Polynomial coefficients (default None) to compute a
902                         gain-elevation correction as a function of
903                         elevation (in degrees).
904            filename:    The name of an ascii file holding correction factors.
905                         The first row of the ascii file must give the column
906                         names and these MUST include columns
907                         "ELEVATION" (degrees) and "FACTOR" (multiply data
908                         by this) somewhere.
909                         The second row must give the data type of the
910                         column. Use 'R' for Real and 'I' for Integer.
911                         An example file would be
912                         (actual factors are arbitrary) :
913
914                         TIME ELEVATION FACTOR
915                         R R R
916                         0.1 0 0.8
917                         0.2 20 0.85
918                         0.3 40 0.9
919                         0.4 60 0.85
920                         0.5 80 0.8
921                         0.6 90 0.75
922            method:      Interpolation method when correcting from a table.
923                         Values are  "nearest", "linear" (default), "cubic"
924                         and "spline"
925            insitu:      if False a new scantable is returned.
926                         Otherwise, the scaling is done in-situ
927                         The default is taken from .asaprc (False)
928            allaxes:     If True apply to all spectra. Otherwise
929                         apply only to the selected (beam/pol/if) spectra only
930                         The default is taken from .asaprc (True if none)
931        """
932
933        if allaxes is None: allaxes = rcParams['scantable.allaxes']
934        if insitu is None: insitu = rcParams['insitu']
935        varlist = vars()
936        if poly is None:
937           poly = ()
938        from os.path import expandvars
939        filename = expandvars(filename)
940        if not insitu:
941            from asap._asap import gainel as _gainEl
942            s = scantable(_gainEl(self, poly, filename, method, allaxes))
943            s._add_history("gain_el", varlist)
944            print_log()
945            return s
946        else:
947            from asap._asap import gainel_insitu as _gainEl
948            _gainEl(self, poly, filename, method, allaxes)
949            self._add_history("gain_el", varlist)
950            print_log()
951            return
952
953    def freq_align(self, reftime=None, method='cubic', perif=False,
954                   insitu=None):
955        """
956        Return a scan where all rows have been aligned in frequency/velocity.
957        The alignment frequency frame (e.g. LSRK) is that set by function
958        set_freqframe.
959        Parameters:
960            reftime:     reference time to align at. By default, the time of
961                         the first row of data is used.
962            method:      Interpolation method for regridding the spectra.
963                         Choose from "nearest", "linear", "cubic" (default)
964                         and "spline"
965            perif:       Generate aligners per freqID (no doppler tracking) or
966                         per IF (scan-based doppler tracking)
967            insitu:      if False a new scantable is returned.
968                         Otherwise, the scaling is done in-situ
969                         The default is taken from .asaprc (False)
970        """
971        if insitu is None: insitu = rcParams['insitu']
972        varlist = vars()
973        if reftime is None: reftime = ''
974        perfreqid = not perif
975        if not insitu:
976            from asap._asap import freq_align as _align
977            s = scantable(_align(self, reftime, method, perfreqid))
978            s._add_history("freq_align", varlist)
979            print_log()
980            return s
981        else:
982            from asap._asap import freq_align_insitu as _align
983            _align(self, reftime, method, perfreqid)
984            self._add_history("freq_align", varlist)
985            print_log()
986            return
987
988    def opacity(self, tau, insitu=None, allaxes=None):
989        """
990        Apply an opacity correction. The data
991        and Tsys are multiplied by the correction factor.
992        Parameters:
993            tau:         Opacity from which the correction factor is
994                         exp(tau*ZD)
995                         where ZD is the zenith-distance
996            insitu:      if False a new scantable is returned.
997                         Otherwise, the scaling is done in-situ
998                         The default is taken from .asaprc (False)
999            allaxes:     if True apply to all spectra. Otherwise
1000                         apply only to the selected (beam/pol/if)spectra only
1001                         The default is taken from .asaprc (True if none)
1002        """
1003        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1004        if insitu is None: insitu = rcParams['insitu']
1005        varlist = vars()
1006        if not insitu:
1007            from asap._asap import opacity as _opacity
1008            s = scantable(_opacity(self, tau, allaxes))
1009            s._add_history("opacity", varlist)
1010            print_log()
1011            return s
1012        else:
1013            from asap._asap import opacity_insitu as _opacity
1014            _opacity(self, tau, allaxes)
1015            self._add_history("opacity", varlist)
1016            print_log()
1017            return
1018
1019    def bin(self, width=5, insitu=None):
1020        """
1021        Return a scan where all spectra have been binned up.
1022            width:       The bin width (default=5) in pixels
1023            insitu:      if False a new scantable is returned.
1024                         Otherwise, the scaling is done in-situ
1025                         The default is taken from .asaprc (False)
1026        """
1027        if insitu is None: insitu = rcParams['insitu']
1028        varlist = vars()
1029        if not insitu:
1030            from asap._asap import bin as _bin
1031            s = scantable(_bin(self, width))
1032            s._add_history("bin",varlist)
1033            print_log()
1034            return s
1035        else:
1036            from asap._asap import bin_insitu as _bin
1037            _bin(self, width)
1038            self._add_history("bin",varlist)
1039            print_log()
1040            return
1041
1042
1043    def resample(self, width=5, method='cubic', insitu=None):
1044        """
1045        Return a scan where all spectra have been binned up
1046            width:       The bin width (default=5) in pixels
1047            method:      Interpolation method when correcting from a table.
1048                         Values are  "nearest", "linear", "cubic" (default)
1049                         and "spline"
1050            insitu:      if False a new scantable is returned.
1051                         Otherwise, the scaling is done in-situ
1052                         The default is taken from .asaprc (False)
1053        """
1054        if insitu is None: insitu = rcParams['insitu']
1055        varlist = vars()
1056        if not insitu:
1057            from asap._asap import resample as _resample
1058            s = scantable(_resample(self, method, width))
1059            s._add_history("resample",varlist)
1060            print_log()
1061            return s
1062        else:
1063            from asap._asap import resample_insitu as _resample
1064            _resample(self, method, width)
1065            self._add_history("resample",varlist)
1066            print_log()
1067            return
1068
1069    def average_pol(self, mask=None, weight='none', insitu=None):
1070        """
1071        Average the Polarisations together.
1072        The polarisation cursor of the output scan is set to 0
1073        Parameters:
1074            mask:        An optional mask defining the region, where the
1075                         averaging will be applied. The output will have all
1076                         specified points masked.
1077            weight:      Weighting scheme. 'none' (default), 'var' (1/var(spec)
1078                         weighted), or 'tsys' (1/Tsys**2 weighted)
1079            insitu:      if False a new scantable is returned.
1080                         Otherwise, the scaling is done in-situ
1081                         The default is taken from .asaprc (False)
1082        """
1083        if insitu is None: insitu = rcParams['insitu']
1084        varlist = vars()
1085
1086        if mask is None:
1087            mask = ()
1088        if not insitu:
1089            from asap._asap import averagepol as _avpol
1090            s = scantable(_avpol(self, mask, weight))
1091            s._add_history("average_pol",varlist)
1092            print_log()
1093            return s
1094        else:
1095            from asap._asap import averagepol_insitu as _avpol
1096            _avpol(self, mask, weight)
1097            self._add_history("average_pol",varlist)
1098            print_log()
1099            return
1100
1101    def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
1102        """
1103        Smooth the spectrum by the specified kernel (conserving flux).
1104        Parameters:
1105            scan:       The input scan
1106            kernel:     The type of smoothing kernel. Select from
1107                        'hanning' (default), 'gaussian' and 'boxcar'.
1108                        The first three characters are sufficient.
1109            width:      The width of the kernel in pixels. For hanning this is
1110                        ignored otherwise it defauls to 5 pixels.
1111                        For 'gaussian' it is the Full Width Half
1112                        Maximum. For 'boxcar' it is the full width.
1113            insitu:     if False a new scantable is returned.
1114                        Otherwise, the scaling is done in-situ
1115                        The default is taken from .asaprc (False)
1116            allaxes:    If True (default) apply to all spectra. Otherwise
1117                        apply only to the selected (beam/pol/if)spectra only
1118                        The default is taken from .asaprc (True if none)
1119        Example:
1120             none
1121        """
1122        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1123        if insitu is None: insitu = rcParams['insitu']
1124        varlist = vars()
1125        if not insitu:
1126            from asap._asap import smooth as _smooth
1127            s = scantable(_smooth(self,kernel,width,allaxes))
1128            s._add_history("smooth", varlist)
1129            print_log()
1130            return s
1131        else:
1132            from asap._asap import smooth_insitu as _smooth
1133            _smooth(self,kernel,width,allaxes)
1134            self._add_history("smooth", varlist)
1135            print_log()
1136            return
1137
1138    def poly_baseline(self, mask=None, order=0, insitu=None, allaxes=None):
1139        """
1140        Return a scan which has been baselined (all rows) by a polynomial.
1141        Parameters:
1142            scan:       a scantable
1143            mask:       an optional mask
1144            order:      the order of the polynomial (default is 0)
1145            insitu:     if False a new scantable is returned.
1146                        Otherwise, the scaling is done in-situ
1147                        The default is taken from .asaprc (False)
1148            allaxes:    If True (default) apply to all spectra. Otherwise
1149                        apply only to the selected (beam/pol/if)spectra only
1150                        The default is taken from .asaprc (True if none)
1151        Example:
1152            # return a scan baselined by a third order polynomial,
1153            # not using a mask
1154            bscan = scan.poly_baseline(order=3)
1155        """
1156        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1157        if insitu is None: insitu = rcParams['insitu']
1158        varlist = vars()
1159        if mask is None:
1160            from numarray import ones
1161            mask = list(ones(self.nchan()))
1162        from asap.asapfitter import fitter
1163        f = fitter()
1164        f.set_scan(self, mask)
1165        f.set_function(poly=order)
1166        sf = f.auto_fit(insitu,allaxes)
1167        if insitu:
1168            self._add_history("poly_baseline", varlist)
1169            print_log()
1170            return
1171        else:
1172            sf._add_history("poly_baseline", varlist)
1173            print_log()
1174            return sf
1175
1176    def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
1177                           threshold=3,insitu=None):
1178        """
1179        Return a scan which has been baselined (all rows) by a polynomial.
1180        Spectral lines are detected first using linefinder and masked out
1181        to avoid them affecting the baseline solution.
1182
1183        Parameters:
1184            mask:       an optional mask retreived from scantable
1185            edge:       an optional number of channel to drop at
1186                        the edge of spectrum. If only one value is
1187                        specified, the same number will be dropped from
1188                        both sides of the spectrum. Default is to keep
1189                        all channels
1190            order:      the order of the polynomial (default is 0)
1191            threshold:  the threshold used by line finder. It is better to
1192                        keep it large as only strong lines affect the
1193                        baseline solution.
1194            insitu:     if False a new scantable is returned.
1195                        Otherwise, the scaling is done in-situ
1196                        The default is taken from .asaprc (False)
1197
1198        Example:
1199            scan2=scan.auto_poly_baseline(order=7)
1200        """
1201        if insitu is None: insitu = rcParams['insitu']
1202        varlist = vars()
1203        from asap.asapfitter import fitter
1204        from asap.asaplinefind import linefinder
1205        from asap import _is_sequence_or_number as _is_valid
1206
1207        if not _is_valid(edge, int):
1208
1209            raise RuntimeError, "Parameter 'edge' has to be an integer or a \
1210            pair of integers specified as a tuple"
1211
1212        # setup fitter
1213        f = fitter()
1214        f.set_function(poly=order)
1215
1216        # setup line finder
1217        fl=linefinder()
1218        fl.set_options(threshold=threshold)
1219
1220        if not insitu:
1221            workscan=self.copy()
1222        else:
1223            workscan=self
1224
1225        sel=workscan.get_cursor()
1226        rows=range(workscan.nrow())
1227        from asap import asaplog
1228        for i in range(workscan.nbeam()):
1229            workscan.setbeam(i)
1230            for j in range(workscan.nif()):
1231                workscan.setif(j)
1232                for k in range(workscan.npol()):
1233                    workscan.setpol(k)
1234                    asaplog.push("Processing:")
1235                    msg = 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
1236                    asaplog.push(msg)
1237                    for iRow in rows:
1238                       fl.set_scan(workscan,mask,edge)
1239                       fl.find_lines(iRow)
1240                       f.set_scan(workscan, fl.get_mask())
1241                       f.x=workscan._getabcissa(iRow)
1242                       f.y=workscan._getspectrum(iRow)
1243                       f.data=None
1244                       f.fit()
1245                       x=f.get_parameters()
1246                       workscan._setspectrum(f.fitter.getresidual(),iRow)
1247        workscan.set_cursor(sel[0],sel[1],sel[2])
1248        if not insitu:
1249            return workscan
1250
1251    def rotate_linpolphase(self, angle, allaxes=None):
1252        """
1253        Rotate the phase of the complex polarization O=Q+iU correlation.
1254        This is always done in situ in the raw data.  So if you call this
1255        function more than once then each call rotates the phase further.
1256        Parameters:
1257            angle:   The angle (degrees) to rotate (add) by.
1258            allaxes: If True apply to all spectra. Otherwise
1259                     apply only to the selected (beam/pol/if)spectra only.
1260                     The default is taken from .asaprc (True if none)
1261        Examples:
1262            scan.rotate_linpolphase(2.3)
1263    """
1264        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1265        varlist = vars()
1266        from asap._asap import _rotate_linpolphase as _rotate
1267        _rotate(self, angle, allaxes)
1268        self._add_history("rotate_linpolphase", varlist)
1269        print_log()
1270        return
1271
1272
1273    def rotate_xyphase(self, angle, allaxes=None):
1274        """
1275        Rotate the phase of the XY correlation.  This is always done in situ
1276        in the data.  So if you call this function more than once
1277        then each call rotates the phase further.
1278        Parameters:
1279            angle:   The angle (degrees) to rotate (add) by.
1280            allaxes: If True apply to all spectra. Otherwise
1281                     apply only to the selected (beam/pol/if)spectra only.
1282                     The default is taken from .asaprc (True if none)
1283        Examples:
1284            scan.rotate_xyphase(2.3)
1285        """
1286        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1287        varlist = vars()
1288        from asap._asap import _rotate_xyphase
1289        _rotate_xyphase(self, angle, allaxes)
1290        self._add_history("rotate_xyphase", varlist)
1291        print_log()
1292        return
1293
1294    def invert_phase(self, allaxes=None):
1295        """
1296        Take the complex conjugate of the complex polarisation cross
1297        correlation.  This is always done in situ in the data.  If
1298        you call this function twice, you end up with the original phase
1299        Parameters:
1300            allaxes: If True apply to all spectra. Otherwise
1301                     apply only to the selected (beam/pol/if)spectra only.
1302                     The default is taken from .asaprc (True if none)
1303        Examples:
1304            scan.invert_phase()
1305        """
1306        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1307        varlist = vars()
1308        from asap._asap import _invert_phase
1309        _invert_phase(self, allaxes)
1310        self._add_history("invert_phase", varlist)
1311        print_log()
1312        return
1313
1314    def swap_ppol(self, allaxes=None):
1315        """
1316        Swaps the order of parallel polarisation products (ie AA, BB
1317        becomes BB, AA) .  This is always done in situ in the data.
1318        If you call this function twice, you end up with the original
1319        order
1320        Parameters:
1321            allaxes: If True apply to all spectra. Otherwise
1322                     apply only to the selected (beam/pol/if)spectra only.
1323                     The default is taken from .asaprc (True if none)
1324        Examples:
1325            scan.swap_pol()
1326        """
1327        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1328        varlist = vars()
1329        from asap._asap import _swap_pol
1330        _swap_pol(self, allaxes)
1331        self._add_history("swap_phase", varlist)
1332        print_log()
1333        return
1334
1335
1336    def add(self, offset, insitu=None, allaxes=None):
1337        """
1338        Return a scan where all spectra have the offset added
1339        Parameters:
1340            offset:      the offset
1341            insitu:      if False a new scantable is returned.
1342                         Otherwise, the scaling is done in-situ
1343                         The default is taken from .asaprc (False)
1344            allaxes:     if True apply to all spectra. Otherwise
1345                         apply only to the selected (beam/pol/if)spectra only
1346                         The default is taken from .asaprc (True if none)
1347        """
1348        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1349        if insitu is None: insitu = rcParams['insitu']
1350        varlist = vars()
1351        if not insitu:
1352            from asap._asap import add as _add
1353            s = scantable(_add(self, offset, allaxes))
1354            s._add_history("add",varlist)
1355            print_log()
1356            return s
1357        else:
1358            from asap._asap import add_insitu as _add
1359            _add(self, offset, allaxes)
1360            self._add_history("add",varlist)
1361            print_log()
1362            return
1363
1364    def scale(self, factor, insitu=None, allaxes=None, tsys=True):
1365        """
1366        Return a scan where all spectra are scaled by the give 'factor'
1367        Parameters:
1368            factor:      the scaling factor
1369            insitu:      if False a new scantable is returned.
1370                         Otherwise, the scaling is done in-situ
1371                         The default is taken from .asaprc (False)
1372            allaxes:     if True apply to all spectra. Otherwise
1373                         apply only to the selected (beam/pol/if)spectra only.
1374                         The default is taken from .asaprc (True if none)
1375            tsys:        if True (default) then apply the operation to Tsys
1376                         as well as the data
1377        """
1378        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1379        if insitu is None: insitu = rcParams['insitu']
1380        varlist = vars()
1381        if not insitu:
1382            from asap._asap import scale as _scale
1383            s = scantable(_scale(self, factor, allaxes, tsys))
1384            s._add_history("scale",varlist)
1385            print_log()
1386            return s
1387        else:
1388            from asap._asap import scale_insitu as _scale
1389            _scale(self, factor, allaxes, tsys)
1390            self._add_history("scale",varlist)
1391            print_log()
1392            return
1393
1394    def auto_quotient(self, mode='suffix', preserve=True):
1395        """
1396        This function allows to build quotients automatically.
1397        It assumes the observation to have the same numer of
1398        "ons" and "offs"
1399        It will support "closest off in time" in the future
1400        Parameters:
1401            mode:           the on/off detection mode; 'suffix' (default)
1402                            'suffix' identifies 'off' scans by the
1403                            trailing '_R' (Mopra/Parkes) or
1404                            '_e'/'_w' (Tid)
1405            preserve:       you can preserve (default) the continuum or
1406                            remove it.  The equations used are
1407                            preserve: Output = Toff * (on/off) - Toff
1408                            remove:   Output = Tref * (on/off) - Ton
1409        """
1410        modes = ["suffix","time"]
1411        if not mode in modes:
1412            print "please provide valid mode. Valid modes are %s" % (modes)
1413            return None
1414        from asap._asap import quotient as _quot
1415        if mode == "suffix":
1416            srcs = self.get_scan("*[^_ewR]")
1417            refs = self.get_scan("*[_ewR]")
1418            if isinstance(srcs,scantable) and isinstance(refs,scantable):
1419                from asap import asaplog
1420                ns,nr = srcs.nrow(),refs.nrow()
1421                msg = "Found %i Off and %i On scans" % (ns,nr)
1422                asaplog.push(msg)
1423                if nr > ns:
1424                    asaplog("Found more Off integrations than On scans - dropping excess Offs.")
1425                    refs = refs.get_scan(range(ns))
1426                print_log()
1427                return scantable(_quot(srcs,refs, preserve))
1428            else:
1429                msg = "Couldn't find any on/off pairs"
1430                if rcParams['verbose']:
1431                    print msg
1432                    return
1433                else:
1434                    raise RuntimeError()
1435        else:
1436            if rcParams['verbose']: print "not yet implemented"
1437            return None
1438
1439    def quotient(self, other, isreference=True, preserve=True):
1440        """
1441        Return the quotient of a 'source' (on) scan and a 'reference' (off)
1442        scan.
1443        The reference can have just one row, even if the signal has many.
1444        Otherwise they must have the same number of rows.
1445        The cursor of the output scan is set to 0
1446        Parameters:
1447            other:          the 'other' scan
1448            isreference:    if the 'other' scan is the reference (default)
1449                            or source
1450            preserve:       you can preserve (default) the continuum or
1451                            remove it.  The equations used are
1452                            preserve: Output = Toff * (on/off) - Toff
1453                            remove:   Output = Tref * (on/off) - Ton
1454        Example:
1455            # src is a scantable for an 'on' scan, ref for an 'off' scantable
1456            q1 = src.quotient(ref)
1457            q2 = ref.quotient(src, isreference=False)
1458            # gives the same result
1459        """
1460        from asap._asap import quotient as _quot
1461        try:
1462            s = None
1463            if isreference:
1464                s = scantable(_quot(self, other, preserve))
1465            else:
1466                s = scantable(_quot(other, self, preserve))
1467            print_log()
1468            return s
1469        except RuntimeError,e:
1470            if rcParams['verbose']:
1471                print e
1472                return
1473            else: raise
1474
1475    def freq_switch(self, insitu=None):
1476        """
1477        Apply frequency switching to the data.
1478        Parameters:
1479            insitu:      if False a new scantable is returned.
1480                         Otherwise, the swictching is done in-situ
1481                         The default is taken from .asaprc (False)
1482        Example:
1483            none
1484        """
1485        if insitu is None: insitu = rcParams['insitu']
1486        varlist = vars()
1487        try:
1488            if insitu:
1489                from asap._asap import _frequency_switch_insitu
1490                _frequency_switch_insitu(self)
1491                self._add_history("freq_switch", varlist)
1492                print_log()
1493                return
1494            else:
1495                from asap._asap import _frequency_switch
1496                sf = scantable(_frequency_switch(self))
1497                sf._add_history("freq_switch", varlist)
1498                print_log()
1499                return sf
1500        except RuntimeError,e:
1501            if rcParams['verbose']: print e
1502            else: raise
1503
1504    def recalc_azel(self):
1505        """
1506        Recalculate the azimuth and elevation for each position.
1507        Parameters:
1508            none
1509        Example:
1510        """
1511        varlist = vars()
1512        self._recalc_azel()
1513        self._add_history("recalc_azel", varlist)
1514        print_log()
1515        return
1516
1517    def __add__(self, other):
1518        varlist = vars()
1519        s = None
1520        if isinstance(other, scantable):
1521            from asap._asap import b_operate as _bop
1522            s = scantable(_bop(self, other, 'add', True))
1523        elif isinstance(other, float):
1524            from asap._asap import add as _add
1525            s = scantable(_add(self, other, True))
1526        else:
1527            raise TypeError("Other input is not a scantable or float value")
1528        s._add_history("operator +", varlist)
1529        print_log()
1530        return s
1531
1532    def __sub__(self, other):
1533        """
1534        implicit on all axes and on Tsys
1535        """
1536        varlist = vars()
1537        s = None
1538        if isinstance(other, scantable):
1539            from asap._asap import b_operate as _bop
1540            s = scantable(_bop(self, other, 'sub', True))
1541        elif isinstance(other, float):
1542            from asap._asap import add as _add
1543            s = scantable(_add(self, -other, True))
1544        else:
1545            raise TypeError("Other input is not a scantable or float value")
1546        s._add_history("operator -", varlist)
1547        print_log()
1548        return s
1549
1550    def __mul__(self, other):
1551        """
1552        implicit on all axes and on Tsys
1553        """
1554        varlist = vars()
1555        s = None
1556        if isinstance(other, scantable):
1557            from asap._asap import b_operate as _bop
1558            s = scantable(_bop(self, other, 'mul', True))
1559        elif isinstance(other, float):
1560            if other == 0.0:
1561                raise ZeroDivisionError("Dividing by zero is not recommended")
1562            from asap._asap import scale as _sca
1563            s = scantable(_sca(self, other, True))
1564        else:
1565            raise TypeError("Other input is not a scantable or float value")
1566        s._add_history("operator *", varlist)
1567        print_log()
1568        return s
1569
1570
1571    def __div__(self, other):
1572        """
1573        implicit on all axes and on Tsys
1574        """
1575        varlist = vars()
1576        s = None
1577        if isinstance(other, scantable):
1578            from asap._asap import b_operate as _bop
1579            s = scantable(_bop(self, other, 'div', True))
1580        elif isinstance(other, float):
1581            if other == 0.0:
1582                raise ZeroDivisionError("Dividing by zero is not recommended")
1583            from asap._asap import scale as _sca
1584            s = scantable(_sca(self, 1.0/other, True))
1585        else:
1586            raise TypeError("Other input is not a scantable or float value")
1587        s._add_history("operator /", varlist)
1588        print_log()
1589        return s
1590
1591    def get_fit(self, row=0):
1592        """
1593        Print or return the stored fits for a row in the scantable
1594        Parameters:
1595            row:    the row which the fit has been applied to.
1596        """
1597        if row > self.nrow():
1598            return
1599        from asap import asapfit
1600        fit = asapfit(self._getfit(row))
1601        if rcParams['verbose']:
1602            print fit
1603            return
1604        else:
1605            return fit.as_dict()
1606
1607    def _add_history(self, funcname, parameters):
1608        # create date
1609        sep = "##"
1610        from datetime import datetime
1611        dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1612        hist = dstr+sep
1613        hist += funcname+sep#cdate+sep
1614        if parameters.has_key('self'): del parameters['self']
1615        for k,v in parameters.iteritems():
1616            if type(v) is dict:
1617                for k2,v2 in v.iteritems():
1618                    hist += k2
1619                    hist += "="
1620                    if isinstance(v2,scantable):
1621                        hist += 'scantable'
1622                    elif k2 == 'mask':
1623                        if isinstance(v2,list) or isinstance(v2,tuple):
1624                            hist += str(self._zip_mask(v2))
1625                        else:
1626                            hist += str(v2)
1627                    else:
1628                        hist += str(v2)
1629            else:
1630                hist += k
1631                hist += "="
1632                if isinstance(v,scantable):
1633                    hist += 'scantable'
1634                elif k == 'mask':
1635                    if isinstance(v,list) or isinstance(v,tuple):
1636                        hist += str(self._zip_mask(v))
1637                    else:
1638                        hist += str(v)
1639                else:
1640                    hist += str(v)
1641            hist += sep
1642        hist = hist[:-2] # remove trailing '##'
1643        self._addhistory(hist)
1644
1645
1646    def _zip_mask(self, mask):
1647        mask = list(mask)
1648        i = 0
1649        segments = []
1650        while mask[i:].count(1):
1651            i += mask[i:].index(1)
1652            if mask[i:].count(0):
1653                j = i + mask[i:].index(0)
1654            else:
1655                j = len(mask)
1656            segments.append([i,j])
1657            i = j
1658        return segments
1659
1660    def _get_ordinate_label(self):
1661        fu = "("+self.get_fluxunit()+")"
1662        import re
1663        lbl = "Intensity"
1664        if re.match(".K.",fu):
1665            lbl = "Brightness Temperature "+ fu
1666        elif re.match(".Jy.",fu):
1667            lbl = "Flux density "+ fu
1668        return lbl
1669
Note: See TracBrowser for help on using the repository browser.