source: trunk/python/scantable.py @ 2429

Last change on this file since 2429 was 2429, checked in by Kana Sugimoto, 12 years ago

New Development: No

JIRA Issue: Yes (CAS-2796)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): asap.scantable

Description: a minor modification in scantable.freq_align.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 139.4 KB
Line 
1"""This module defines the scantable class."""
2
3import os
4import tempfile
5import numpy
6try:
7    from functools import wraps as wraps_dec
8except ImportError:
9    from asap.compatibility import wraps as wraps_dec
10
11from asap.env import is_casapy
12from asap._asap import Scantable
13from asap._asap import filler, msfiller
14from asap.parameters import rcParams
15from asap.logging import asaplog, asaplog_post_dec
16from asap.selector import selector
17from asap.linecatalog import linecatalog
18from asap.coordinate import coordinate
19from asap.utils import _n_bools, mask_not, mask_and, mask_or, page
20from asap.asapfitter import fitter
21
22###############################################################
23### WK temporarily added these lines for testing 2011/11/28 ###
24###############################################################
25#from asap._asap import TestClass
26
27#class testclass(TestClass):
28#    def __init__(self, nelem):
29#        TestClass.__init__(self, nelem)
30
31###############################################################
32
33
34def preserve_selection(func):
35    @wraps_dec(func)
36    def wrap(obj, *args, **kw):
37        basesel = obj.get_selection()
38        try:
39            val = func(obj, *args, **kw)
40        finally:
41            obj.set_selection(basesel)
42        return val
43    return wrap
44
45def is_scantable(filename):
46    """Is the given file a scantable?
47
48    Parameters:
49
50        filename: the name of the file/directory to test
51
52    """
53    if ( os.path.isdir(filename)
54         and os.path.exists(filename+'/table.info')
55         and os.path.exists(filename+'/table.dat') ):
56        f=open(filename+'/table.info')
57        l=f.readline()
58        f.close()
59        #if ( l.find('Scantable') != -1 ):
60        if ( l.find('Measurement Set') == -1 ):
61            return True
62        else:
63            return False
64    else:
65        return False
66##     return (os.path.isdir(filename)
67##             and not os.path.exists(filename+'/table.f1')
68##             and os.path.exists(filename+'/table.info'))
69
70def is_ms(filename):
71    """Is the given file a MeasurementSet?
72
73    Parameters:
74
75        filename: the name of the file/directory to test
76
77    """
78    if ( os.path.isdir(filename)
79         and os.path.exists(filename+'/table.info')
80         and os.path.exists(filename+'/table.dat') ):
81        f=open(filename+'/table.info')
82        l=f.readline()
83        f.close()
84        if ( l.find('Measurement Set') != -1 ):
85            return True
86        else:
87            return False
88    else:
89        return False
90
91def normalise_edge_param(edge):
92    """\
93    Convert a given edge value to a one-dimensional array that can be
94    given to baseline-fitting/subtraction functions.
95    The length of the output value will be an even because values for
96    the both sides of spectra are to be contained for each IF. When
97    the length is 2, the values will be applied to all IFs. If the length
98    is larger than 2, it will be 2*ifnos().
99    Accepted format of edge include:
100            * an integer - will be used for both sides of spectra of all IFs.
101                  e.g. 10 is converted to [10,10]
102            * an empty list/tuple [] - converted to [0, 0] and used for all IFs.
103            * a list/tuple containing an integer - same as the above case.
104                  e.g. [10] is converted to [10,10]
105            * a list/tuple containing two integers - will be used for all IFs.
106                  e.g. [5,10] is output as it is. no need to convert.
107            * a list/tuple of lists/tuples containing TWO integers -
108                  each element of edge will be used for each IF.
109                  e.g. [[5,10],[15,20]] - [5,10] for IF[0] and [15,20] for IF[1].
110                 
111                  If an element contains the same integer values, the input 'edge'
112                  parameter can be given in a simpler shape in the following cases:
113                      ** when len(edge)!=2
114                          any elements containing the same values can be replaced
115                          to single integers.
116                          e.g. [[15,15]] can be simplified to [15] (or [15,15] or 15 also).
117                          e.g. [[1,1],[2,2],[3,3]] can be simplified to [1,2,3].
118                      ** when len(edge)=2
119                          care is needed for this case: ONLY ONE of the
120                          elements can be a single integer,
121                          e.g. [[5,5],[10,10]] can be simplified to [5,[10,10]]
122                               or [[5,5],10], but can NOT be simplified to [5,10].
123                               when [5,10] given, it is interpreted as
124                               [[5,10],[5,10],[5,10],...] instead, as shown before.
125    """
126    from asap import _is_sequence_or_number as _is_valid
127    if isinstance(edge, list) or isinstance(edge, tuple):
128        for edgepar in edge:
129            if not _is_valid(edgepar, int):
130                raise ValueError, "Each element of the 'edge' tuple has \
131                                   to be a pair of integers or an integer."
132            if isinstance(edgepar, list) or isinstance(edgepar, tuple):
133                if len(edgepar) != 2:
134                    raise ValueError, "Each element of the 'edge' tuple has \
135                                       to be a pair of integers or an integer."
136    else:
137        if not _is_valid(edge, int):
138            raise ValueError, "Parameter 'edge' has to be an integer or a \
139                               pair of integers specified as a tuple. \
140                               Nested tuples are allowed \
141                               to make individual selection for different IFs."
142       
143
144    if isinstance(edge, int):
145        edge = [ edge, edge ]                 # e.g. 3   => [3,3]
146    elif isinstance(edge, list) or isinstance(edge, tuple):
147        if len(edge) == 0:
148            edge = [0, 0]                     # e.g. []  => [0,0]
149        elif len(edge) == 1:
150            if isinstance(edge[0], int):
151                edge = [ edge[0], edge[0] ]   # e.g. [1] => [1,1]
152
153    commonedge = True
154    if len(edge) > 2: commonedge = False
155    else:
156        for edgepar in edge:
157            if isinstance(edgepar, list) or isinstance(edgepar, tuple):
158                commonedge = False
159                break
160
161    if commonedge:
162        if len(edge) > 1:
163            norm_edge = edge
164        else:
165            norm_edge = edge + edge
166    else:
167        norm_edge = []
168        for edgepar in edge:
169            if isinstance(edgepar, int):
170                norm_edge += [edgepar, edgepar]
171            else:
172                norm_edge += edgepar
173
174    return norm_edge
175
176def raise_fitting_failure_exception(e):
177    msg = "The fit failed, possibly because it didn't converge."
178    if rcParams["verbose"]:
179        asaplog.push(str(e))
180        asaplog.push(str(msg))
181    else:
182        raise RuntimeError(str(e)+'\n'+msg)
183
184def pack_progress_params(showprogress, minnrow):
185    return str(showprogress).lower() + ',' + str(minnrow)
186
187class scantable(Scantable):
188    """\
189        The ASAP container for scans (single-dish data).
190    """
191
192    @asaplog_post_dec
193    def __init__(self, filename, average=None, unit=None, parallactify=None,
194                 **args):
195        """\
196        Create a scantable from a saved one or make a reference
197
198        Parameters:
199
200            filename:     the name of an asap table on disk
201                          or
202                          the name of a rpfits/sdfits/ms file
203                          (integrations within scans are auto averaged
204                          and the whole file is read) or
205                          [advanced] a reference to an existing scantable
206
207            average:      average all integrations withinb a scan on read.
208                          The default (True) is taken from .asaprc.
209
210            unit:         brightness unit; must be consistent with K or Jy.
211                          Over-rides the default selected by the filler
212                          (input rpfits/sdfits/ms) or replaces the value
213                          in existing scantables
214
215            antenna:      for MeasurementSet input data only:
216                          Antenna selection. integer (id) or string
217                          (name or id).
218
219            parallactify: Indicate that the data had been parallactified.
220                          Default (false) is taken from rc file.
221
222        """
223        if average is None:
224            average = rcParams['scantable.autoaverage']
225        parallactify = parallactify or rcParams['scantable.parallactify']
226        varlist = vars()
227        from asap._asap import stmath
228        self._math = stmath( rcParams['insitu'] )
229        if isinstance(filename, Scantable):
230            Scantable.__init__(self, filename)
231        else:
232            if isinstance(filename, str):
233                filename = os.path.expandvars(filename)
234                filename = os.path.expanduser(filename)
235                if not os.path.exists(filename):
236                    s = "File '%s' not found." % (filename)
237                    raise IOError(s)
238                if is_scantable(filename):
239                    ondisk = rcParams['scantable.storage'] == 'disk'
240                    Scantable.__init__(self, filename, ondisk)
241                    if unit is not None:
242                        self.set_fluxunit(unit)
243                    if average:
244                        self._assign( self.average_time( scanav=True ) )
245                    # do not reset to the default freqframe
246                    #self.set_freqframe(rcParams['scantable.freqframe'])
247                elif is_ms(filename):
248                    # Measurement Set
249                    opts={'ms': {}}
250                    mskeys=['getpt','antenna']
251                    for key in mskeys:
252                        if key in args.keys():
253                            opts['ms'][key] = args[key]
254                    self._fill([filename], unit, average, opts)
255                elif os.path.isfile(filename):
256                    self._fill([filename], unit, average)
257                    # only apply to new data not "copy constructor"
258                    self.parallactify(parallactify)
259                else:
260                    msg = "The given file '%s'is not a valid " \
261                          "asap table." % (filename)
262                    raise IOError(msg)
263            elif (isinstance(filename, list) or isinstance(filename, tuple)) \
264                  and isinstance(filename[-1], str):
265                self._fill(filename, unit, average)
266        self.parallactify(parallactify)
267        self._add_history("scantable", varlist)
268
269    @asaplog_post_dec
270    def save(self, name=None, format=None, overwrite=False):
271        """\
272        Store the scantable on disk. This can be an asap (aips++) Table,
273        SDFITS or MS2 format.
274
275        Parameters:
276
277            name:        the name of the outputfile. For format "ASCII"
278                         this is the root file name (data in 'name'.txt
279                         and header in 'name'_header.txt)
280
281            format:      an optional file format. Default is ASAP.
282                         Allowed are:
283
284                            * 'ASAP' (save as ASAP [aips++] Table),
285                            * 'SDFITS' (save as SDFITS file)
286                            * 'ASCII' (saves as ascii text file)
287                            * 'MS2' (saves as an casacore MeasurementSet V2)
288                            * 'FITS' (save as image FITS - not readable by
289                                      class)
290                            * 'CLASS' (save as FITS readable by CLASS)
291
292            overwrite:   If the file should be overwritten if it exists.
293                         The default False is to return with warning
294                         without writing the output. USE WITH CARE.
295
296        Example::
297
298            scan.save('myscan.asap')
299            scan.save('myscan.sdfits', 'SDFITS')
300
301        """
302        from os import path
303        format = format or rcParams['scantable.save']
304        suffix = '.'+format.lower()
305        if name is None or name == "":
306            name = 'scantable'+suffix
307            msg = "No filename given. Using default name %s..." % name
308            asaplog.push(msg)
309        name = path.expandvars(name)
310        if path.isfile(name) or path.isdir(name):
311            if not overwrite:
312                msg = "File %s exists." % name
313                raise IOError(msg)
314        format2 = format.upper()
315        if format2 == 'ASAP':
316            self._save(name)
317        elif format2 == 'MS2':
318            msopt = {'ms': {'overwrite': overwrite } }
319            from asap._asap import mswriter
320            writer = mswriter( self )
321            writer.write( name, msopt )
322        else:
323            from asap._asap import stwriter as stw
324            writer = stw(format2)
325            writer.write(self, name)
326        return
327
328    def copy(self):
329        """Return a copy of this scantable.
330
331        *Note*:
332
333            This makes a full (deep) copy. scan2 = scan1 makes a reference.
334
335        Example::
336
337            copiedscan = scan.copy()
338
339        """
340        sd = scantable(Scantable._copy(self))
341        return sd
342
343    def drop_scan(self, scanid=None):
344        """\
345        Return a new scantable where the specified scan number(s) has(have)
346        been dropped.
347
348        Parameters:
349
350            scanid:    a (list of) scan number(s)
351
352        """
353        from asap import _is_sequence_or_number as _is_valid
354        from asap import _to_list
355        from asap import unique
356        if not _is_valid(scanid):
357            raise RuntimeError( 'Please specify a scanno to drop from the'
358                                ' scantable' )
359        scanid = _to_list(scanid)
360        allscans = unique([ self.getscan(i) for i in range(self.nrow())])
361        for sid in scanid: allscans.remove(sid)
362        if len(allscans) == 0:
363            raise ValueError("Can't remove all scans")
364        sel = selector(scans=allscans)
365        return self._select_copy(sel)
366
367    def _select_copy(self, selection):
368        orig = self.get_selection()
369        self.set_selection(orig+selection)
370        cp = self.copy()
371        self.set_selection(orig)
372        return cp
373
374    def get_scan(self, scanid=None):
375        """\
376        Return a specific scan (by scanno) or collection of scans (by
377        source name) in a new scantable.
378
379        *Note*:
380
381            See scantable.drop_scan() for the inverse operation.
382
383        Parameters:
384
385            scanid:    a (list of) scanno or a source name, unix-style
386                       patterns are accepted for source name matching, e.g.
387                       '*_R' gets all 'ref scans
388
389        Example::
390
391            # get all scans containing the source '323p459'
392            newscan = scan.get_scan('323p459')
393            # get all 'off' scans
394            refscans = scan.get_scan('*_R')
395            # get a susbset of scans by scanno (as listed in scan.summary())
396            newscan = scan.get_scan([0, 2, 7, 10])
397
398        """
399        if scanid is None:
400            raise RuntimeError( 'Please specify a scan no or name to '
401                                'retrieve from the scantable' )
402        try:
403            bsel = self.get_selection()
404            sel = selector()
405            if type(scanid) is str:
406                sel.set_name(scanid)
407                return self._select_copy(sel)
408            elif type(scanid) is int:
409                sel.set_scans([scanid])
410                return self._select_copy(sel)
411            elif type(scanid) is list:
412                sel.set_scans(scanid)
413                return self._select_copy(sel)
414            else:
415                msg = "Illegal scanid type, use 'int' or 'list' if ints."
416                raise TypeError(msg)
417        except RuntimeError:
418            raise
419
420    def __str__(self):
421        tempFile = tempfile.NamedTemporaryFile()
422        Scantable._summary(self, tempFile.name)
423        tempFile.seek(0)
424        asaplog.clear()
425        return tempFile.file.read()
426
427    @asaplog_post_dec
428    def summary(self, filename=None):
429        """\
430        Print a summary of the contents of this scantable.
431
432        Parameters:
433
434            filename:    the name of a file to write the putput to
435                         Default - no file output
436
437        """
438#         info = Scantable._summary(self)
439        if filename is not None:
440            if filename is "":
441                filename = 'scantable_summary.txt'
442            from os.path import expandvars, isdir
443            filename = expandvars(filename)
444#             if not isdir(filename):
445#                 data = open(filename, 'w')
446#                 data.write(info)
447#                 data.close()
448#             else:
449            if isdir(filename):
450                msg = "Illegal file name '%s'." % (filename)
451                raise IOError(msg)
452        else:
453            filename = ""
454        Scantable._summary(self, filename)
455#         info = Scantable._summary(self, filename)
456#         return page(info)
457
458    def get_spectrum(self, rowno):
459        """Return the spectrum for the current row in the scantable as a list.
460
461        Parameters:
462
463             rowno:   the row number to retrieve the spectrum from
464
465        """
466        return self._getspectrum(rowno)
467
468    def get_mask(self, rowno):
469        """Return the mask for the current row in the scantable as a list.
470
471        Parameters:
472
473             rowno:   the row number to retrieve the mask from
474
475        """
476        return self._getmask(rowno)
477
478    def set_spectrum(self, spec, rowno):
479        """Set the spectrum for the current row in the scantable.
480
481        Parameters:
482
483             spec:   the new spectrum
484
485             rowno:  the row number to set the spectrum for
486
487        """
488        assert(len(spec) == self.nchan(self.getif(rowno)))
489        return self._setspectrum(spec, rowno)
490
491    def get_coordinate(self, rowno):
492        """Return the (spectral) coordinate for a a given 'rowno'.
493
494        *Note*:
495
496            * This coordinate is only valid until a scantable method modifies
497              the frequency axis.
498            * This coordinate does contain the original frequency set-up
499              NOT the new frame. The conversions however are done using the user
500              specified frame (e.g. LSRK/TOPO). To get the 'real' coordinate,
501              use scantable.freq_align first. Without it there is no closure,
502              i.e.::
503
504                  c = myscan.get_coordinate(0)
505                  c.to_frequency(c.get_reference_pixel()) != c.get_reference_value()
506
507        Parameters:
508
509             rowno:    the row number for the spectral coordinate
510
511        """
512        return coordinate(Scantable.get_coordinate(self, rowno))
513
514    def get_selection(self):
515        """\
516        Get the selection object currently set on this scantable.
517
518        Example::
519
520            sel = scan.get_selection()
521            sel.set_ifs(0)              # select IF 0
522            scan.set_selection(sel)     # apply modified selection
523
524        """
525        return selector(self._getselection())
526
527    def set_selection(self, selection=None, **kw):
528        """\
529        Select a subset of the data. All following operations on this scantable
530        are only applied to thi selection.
531
532        Parameters:
533
534            selection:    a selector object (default unset the selection), or
535                          any combination of "pols", "ifs", "beams", "scans",
536                          "cycles", "name", "query"
537
538        Examples::
539
540            sel = selector()         # create a selection object
541            self.set_scans([0, 3])    # select SCANNO 0 and 3
542            scan.set_selection(sel)  # set the selection
543            scan.summary()           # will only print summary of scanno 0 an 3
544            scan.set_selection()     # unset the selection
545            # or the equivalent
546            scan.set_selection(scans=[0,3])
547            scan.summary()           # will only print summary of scanno 0 an 3
548            scan.set_selection()     # unset the selection
549
550        """
551        if selection is None:
552            # reset
553            if len(kw) == 0:
554                selection = selector()
555            else:
556                # try keywords
557                for k in kw:
558                    if k not in selector.fields:
559                        raise KeyError("Invalid selection key '%s', "
560                                       "valid keys are %s" % (k,
561                                                              selector.fields))
562                selection = selector(**kw)
563        self._setselection(selection)
564
565    def get_row(self, row=0, insitu=None):
566        """\
567        Select a row in the scantable.
568        Return a scantable with single row.
569
570        Parameters:
571
572            row:    row no of integration, default is 0.
573            insitu: if False a new scantable is returned. Otherwise, the
574                    scaling is done in-situ. The default is taken from .asaprc
575                    (False)
576
577        """
578        if insitu is None:
579            insitu = rcParams['insitu']
580        if not insitu:
581            workscan = self.copy()
582        else:
583            workscan = self
584        # Select a row
585        sel = selector()
586        sel.set_rows([row])
587        workscan.set_selection(sel)
588        if not workscan.nrow() == 1:
589            msg = "Could not identify single row. %d rows selected." \
590                % (workscan.nrow())
591            raise RuntimeError(msg)
592        if insitu:
593            self._assign(workscan)
594        else:
595            return workscan
596
597    @asaplog_post_dec
598    def stats(self, stat='stddev', mask=None, form='3.3f', row=None):
599        """\
600        Determine the specified statistic of the current beam/if/pol
601        Takes a 'mask' as an optional parameter to specify which
602        channels should be excluded.
603
604        Parameters:
605
606            stat:    'min', 'max', 'min_abc', 'max_abc', 'sumsq', 'sum',
607                     'mean', 'var', 'stddev', 'avdev', 'rms', 'median'
608
609            mask:    an optional mask specifying where the statistic
610                     should be determined.
611
612            form:    format string to print statistic values
613
614            row:     row number of spectrum to process.
615                     (default is None: for all rows)
616
617        Example:
618            scan.set_unit('channel')
619            msk = scan.create_mask([100, 200], [500, 600])
620            scan.stats(stat='mean', mask=m)
621
622        """
623        mask = mask or []
624        if not self._check_ifs():
625            raise ValueError("Cannot apply mask as the IFs have different "
626                             "number of channels. Please use setselection() "
627                             "to select individual IFs")
628        rtnabc = False
629        if stat.lower().endswith('_abc'): rtnabc = True
630        getchan = False
631        if stat.lower().startswith('min') or stat.lower().startswith('max'):
632            chan = self._math._minmaxchan(self, mask, stat)
633            getchan = True
634            statvals = []
635        if not rtnabc:
636            if row == None:
637                statvals = self._math._stats(self, mask, stat)
638            else:
639                statvals = self._math._statsrow(self, mask, stat, int(row))
640
641        #def cb(i):
642        #    return statvals[i]
643
644        #return self._row_callback(cb, stat)
645
646        label=stat
647        #callback=cb
648        out = ""
649        #outvec = []
650        sep = '-'*50
651
652        if row == None:
653            rows = xrange(self.nrow())
654        elif isinstance(row, int):
655            rows = [ row ]
656
657        for i in rows:
658            refstr = ''
659            statunit= ''
660            if getchan:
661                qx, qy = self.chan2data(rowno=i, chan=chan[i])
662                if rtnabc:
663                    statvals.append(qx['value'])
664                    refstr = ('(value: %'+form) % (qy['value'])+' ['+qy['unit']+'])'
665                    statunit= '['+qx['unit']+']'
666                else:
667                    refstr = ('(@ %'+form) % (qx['value'])+' ['+qx['unit']+'])'
668
669            tm = self._gettime(i)
670            src = self._getsourcename(i)
671            out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
672            out += 'Time[%s]:\n' % (tm)
673            if self.nbeam(-1) > 1: out +=  ' Beam[%d] ' % (self.getbeam(i))
674            if self.nif(-1) > 1:   out +=  ' IF[%d] ' % (self.getif(i))
675            if self.npol(-1) > 1:  out +=  ' Pol[%d] ' % (self.getpol(i))
676            #outvec.append(callback(i))
677            if len(rows) > 1:
678                # out += ('= %'+form) % (outvec[i]) +'   '+refstr+'\n'
679                out += ('= %'+form) % (statvals[i]) +'   '+refstr+'\n'
680            else:
681                # out += ('= %'+form) % (outvec[0]) +'   '+refstr+'\n'
682                out += ('= %'+form) % (statvals[0]) +'   '+refstr+'\n'
683            out +=  sep+"\n"
684
685        import os
686        if os.environ.has_key( 'USER' ):
687            usr = os.environ['USER']
688        else:
689            import commands
690            usr = commands.getoutput( 'whoami' )
691        tmpfile = '/tmp/tmp_'+usr+'_casapy_asap_scantable_stats'
692        f = open(tmpfile,'w')
693        print >> f, sep
694        print >> f, ' %s %s' % (label, statunit)
695        print >> f, sep
696        print >> f, out
697        f.close()
698        f = open(tmpfile,'r')
699        x = f.readlines()
700        f.close()
701        asaplog.push(''.join(x), False)
702
703        return statvals
704
705    def chan2data(self, rowno=0, chan=0):
706        """\
707        Returns channel/frequency/velocity and spectral value
708        at an arbitrary row and channel in the scantable.
709
710        Parameters:
711
712            rowno:   a row number in the scantable. Default is the
713                     first row, i.e. rowno=0
714
715            chan:    a channel in the scantable. Default is the first
716                     channel, i.e. pos=0
717
718        """
719        if isinstance(rowno, int) and isinstance(chan, int):
720            qx = {'unit': self.get_unit(),
721                  'value': self._getabcissa(rowno)[chan]}
722            qy = {'unit': self.get_fluxunit(),
723                  'value': self._getspectrum(rowno)[chan]}
724            return qx, qy
725
726    def stddev(self, mask=None):
727        """\
728        Determine the standard deviation of the current beam/if/pol
729        Takes a 'mask' as an optional parameter to specify which
730        channels should be excluded.
731
732        Parameters:
733
734            mask:    an optional mask specifying where the standard
735                     deviation should be determined.
736
737        Example::
738
739            scan.set_unit('channel')
740            msk = scan.create_mask([100, 200], [500, 600])
741            scan.stddev(mask=m)
742
743        """
744        return self.stats(stat='stddev', mask=mask);
745
746
747    def get_column_names(self):
748        """\
749        Return a  list of column names, which can be used for selection.
750        """
751        return list(Scantable.get_column_names(self))
752
753    def get_tsys(self, row=-1):
754        """\
755        Return the System temperatures.
756
757        Parameters:
758
759            row:    the rowno to get the information for. (default all rows)
760
761        Returns:
762
763            a list of Tsys values for the current selection
764
765        """
766        if row > -1:
767            return self._get_column(self._gettsys, row)
768        return self._row_callback(self._gettsys, "Tsys")
769
770    def get_tsysspectrum(self, row=-1):
771        """\
772        Return the channel dependent system temperatures.
773
774        Parameters:
775
776            row:    the rowno to get the information for. (default all rows)
777
778        Returns:
779
780            a list of Tsys values for the current selection
781
782        """
783        return self._get_column( self._gettsysspectrum, row )
784
785    def get_weather(self, row=-1):
786        """\
787        Return the weather informations.
788
789        Parameters:
790
791            row:    the rowno to get the information for. (default all rows)
792
793        Returns:
794
795            a dict or list of of dicts of values for the current selection
796
797        """
798
799        values = self._get_column(self._get_weather, row)
800        if row > -1:
801            return {'temperature': values[0],
802                    'pressure': values[1], 'humidity' : values[2],
803                    'windspeed' : values[3], 'windaz' : values[4]
804                    }
805        else:
806            out = []
807            for r in values:
808
809                out.append({'temperature': r[0],
810                            'pressure': r[1], 'humidity' : r[2],
811                            'windspeed' : r[3], 'windaz' : r[4]
812                    })
813            return out
814
815    def _row_callback(self, callback, label):
816        out = ""
817        outvec = []
818        sep = '-'*50
819        for i in range(self.nrow()):
820            tm = self._gettime(i)
821            src = self._getsourcename(i)
822            out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
823            out += 'Time[%s]:\n' % (tm)
824            if self.nbeam(-1) > 1:
825                out +=  ' Beam[%d] ' % (self.getbeam(i))
826            if self.nif(-1) > 1: out +=  ' IF[%d] ' % (self.getif(i))
827            if self.npol(-1) > 1: out +=  ' Pol[%d] ' % (self.getpol(i))
828            outvec.append(callback(i))
829            out += '= %3.3f\n' % (outvec[i])
830            out +=  sep+'\n'
831
832        asaplog.push(sep)
833        asaplog.push(" %s" % (label))
834        asaplog.push(sep)
835        asaplog.push(out)
836        asaplog.post()
837        return outvec
838
839    def _get_column(self, callback, row=-1, *args):
840        """
841        """
842        if row == -1:
843            return [callback(i, *args) for i in range(self.nrow())]
844        else:
845            if  0 <= row < self.nrow():
846                return callback(row, *args)
847
848
849    def get_time(self, row=-1, asdatetime=False, prec=-1):
850        """\
851        Get a list of time stamps for the observations.
852        Return a datetime object or a string (default) for each
853        integration time stamp in the scantable.
854
855        Parameters:
856
857            row:          row no of integration. Default -1 return all rows
858
859            asdatetime:   return values as datetime objects rather than strings
860
861            prec:         number of digits shown. Default -1 to automatic
862                          calculation.
863                          Note this number is equals to the digits of MVTime,
864                          i.e., 0<prec<3: dates with hh:: only,
865                          <5: with hh:mm:, <7 or 0: with hh:mm:ss,
866                          and 6> : with hh:mm:ss.tt... (prec-6 t's added)
867
868        """
869        from datetime import datetime
870        if prec < 0:
871            # automagically set necessary precision +1
872            prec = 7 - \
873                numpy.floor(numpy.log10(numpy.min(self.get_inttime(row))))
874            prec = max(6, int(prec))
875        else:
876            prec = max(0, prec)
877        if asdatetime:
878            #precision can be 1 millisecond at max
879            prec = min(12, prec)
880
881        times = self._get_column(self._gettime, row, prec)
882        if not asdatetime:
883            return times
884        format = "%Y/%m/%d/%H:%M:%S.%f"
885        if prec < 7:
886            nsub = 1 + (((6-prec)/2) % 3)
887            substr = [".%f","%S","%M"]
888            for i in range(nsub):
889                format = format.replace(substr[i],"")
890        if isinstance(times, list):
891            return [datetime.strptime(i, format) for i in times]
892        else:
893            return datetime.strptime(times, format)
894
895
896    def get_inttime(self, row=-1):
897        """\
898        Get a list of integration times for the observations.
899        Return a time in seconds for each integration in the scantable.
900
901        Parameters:
902
903            row:    row no of integration. Default -1 return all rows.
904
905        """
906        return self._get_column(self._getinttime, row)
907
908
909    def get_sourcename(self, row=-1):
910        """\
911        Get a list source names for the observations.
912        Return a string for each integration in the scantable.
913        Parameters:
914
915            row:    row no of integration. Default -1 return all rows.
916
917        """
918        return self._get_column(self._getsourcename, row)
919
920    def get_elevation(self, row=-1):
921        """\
922        Get a list of elevations for the observations.
923        Return a float for each integration in the scantable.
924
925        Parameters:
926
927            row:    row no of integration. Default -1 return all rows.
928
929        """
930        return self._get_column(self._getelevation, row)
931
932    def get_azimuth(self, row=-1):
933        """\
934        Get a list of azimuths for the observations.
935        Return a float for each integration in the scantable.
936
937        Parameters:
938            row:    row no of integration. Default -1 return all rows.
939
940        """
941        return self._get_column(self._getazimuth, row)
942
943    def get_parangle(self, row=-1):
944        """\
945        Get a list of parallactic angles for the observations.
946        Return a float for each integration in the scantable.
947
948        Parameters:
949
950            row:    row no of integration. Default -1 return all rows.
951
952        """
953        return self._get_column(self._getparangle, row)
954
955    def get_direction(self, row=-1):
956        """
957        Get a list of Positions on the sky (direction) for the observations.
958        Return a string for each integration in the scantable.
959
960        Parameters:
961
962            row:    row no of integration. Default -1 return all rows
963
964        """
965        return self._get_column(self._getdirection, row)
966
967    def get_directionval(self, row=-1):
968        """\
969        Get a list of Positions on the sky (direction) for the observations.
970        Return a float for each integration in the scantable.
971
972        Parameters:
973
974            row:    row no of integration. Default -1 return all rows
975
976        """
977        return self._get_column(self._getdirectionvec, row)
978
979    @asaplog_post_dec
980    def set_unit(self, unit='channel'):
981        """\
982        Set the unit for all following operations on this scantable
983
984        Parameters:
985
986            unit:    optional unit, default is 'channel'. Use one of '*Hz',
987                     'km/s', 'channel' or equivalent ''
988
989        """
990        varlist = vars()
991        if unit in ['', 'pixel', 'channel']:
992            unit = ''
993        inf = list(self._getcoordinfo())
994        inf[0] = unit
995        self._setcoordinfo(inf)
996        self._add_history("set_unit", varlist)
997
998    @asaplog_post_dec
999    def set_instrument(self, instr):
1000        """\
1001        Set the instrument for subsequent processing.
1002
1003        Parameters:
1004
1005            instr:    Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
1006                      'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
1007
1008        """
1009        self._setInstrument(instr)
1010        self._add_history("set_instument", vars())
1011
1012    @asaplog_post_dec
1013    def set_feedtype(self, feedtype):
1014        """\
1015        Overwrite the feed type, which might not be set correctly.
1016
1017        Parameters:
1018
1019            feedtype:     'linear' or 'circular'
1020
1021        """
1022        self._setfeedtype(feedtype)
1023        self._add_history("set_feedtype", vars())
1024
1025    @asaplog_post_dec
1026    def set_doppler(self, doppler='RADIO'):
1027        """\
1028        Set the doppler for all following operations on this scantable.
1029
1030        Parameters:
1031
1032            doppler:    One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
1033
1034        """
1035        varlist = vars()
1036        inf = list(self._getcoordinfo())
1037        inf[2] = doppler
1038        self._setcoordinfo(inf)
1039        self._add_history("set_doppler", vars())
1040
1041    @asaplog_post_dec
1042    def set_freqframe(self, frame=None):
1043        """\
1044        Set the frame type of the Spectral Axis.
1045
1046        Parameters:
1047
1048            frame:   an optional frame type, default 'LSRK'. Valid frames are:
1049                     'TOPO', 'LSRD', 'LSRK', 'BARY',
1050                     'GEO', 'GALACTO', 'LGROUP', 'CMB'
1051
1052        Example::
1053
1054            scan.set_freqframe('BARY')
1055
1056        """
1057        frame = frame or rcParams['scantable.freqframe']
1058        varlist = vars()
1059        # "REST" is not implemented in casacore
1060        #valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \
1061        #           'GEO', 'GALACTO', 'LGROUP', 'CMB']
1062        valid = ['TOPO', 'LSRD', 'LSRK', 'BARY', \
1063                   'GEO', 'GALACTO', 'LGROUP', 'CMB']
1064
1065        if frame in valid:
1066            inf = list(self._getcoordinfo())
1067            inf[1] = frame
1068            self._setcoordinfo(inf)
1069            self._add_history("set_freqframe", varlist)
1070        else:
1071            msg  = "Please specify a valid freq type. Valid types are:\n", valid
1072            raise TypeError(msg)
1073
1074    @asaplog_post_dec
1075    def set_dirframe(self, frame=""):
1076        """\
1077        Set the frame type of the Direction on the sky.
1078
1079        Parameters:
1080
1081            frame:   an optional frame type, default ''. Valid frames are:
1082                     'J2000', 'B1950', 'GALACTIC'
1083
1084        Example:
1085
1086            scan.set_dirframe('GALACTIC')
1087
1088        """
1089        varlist = vars()
1090        Scantable.set_dirframe(self, frame)
1091        self._add_history("set_dirframe", varlist)
1092
1093    def get_unit(self):
1094        """\
1095        Get the default unit set in this scantable
1096
1097        Returns:
1098
1099            A unit string
1100
1101        """
1102        inf = self._getcoordinfo()
1103        unit = inf[0]
1104        if unit == '': unit = 'channel'
1105        return unit
1106
1107    @asaplog_post_dec
1108    def get_abcissa(self, rowno=0):
1109        """\
1110        Get the abcissa in the current coordinate setup for the currently
1111        selected Beam/IF/Pol
1112
1113        Parameters:
1114
1115            rowno:    an optional row number in the scantable. Default is the
1116                      first row, i.e. rowno=0
1117
1118        Returns:
1119
1120            The abcissa values and the format string (as a dictionary)
1121
1122        """
1123        abc = self._getabcissa(rowno)
1124        lbl = self._getabcissalabel(rowno)
1125        return abc, lbl
1126
1127    @asaplog_post_dec
1128    def flag(self, mask=None, unflag=False, row=-1):
1129        """\
1130        Flag the selected data using an optional channel mask.
1131
1132        Parameters:
1133
1134            mask:   an optional channel mask, created with create_mask. Default
1135                    (no mask) is all channels.
1136
1137            unflag:    if True, unflag the data
1138
1139            row:    an optional row number in the scantable.
1140                      Default -1 flags all rows
1141                     
1142        """
1143        varlist = vars()
1144        mask = mask or []
1145        self._flag(row, mask, unflag)
1146        self._add_history("flag", varlist)
1147
1148    @asaplog_post_dec
1149    def flag_row(self, rows=None, unflag=False):
1150        """\
1151        Flag the selected data in row-based manner.
1152
1153        Parameters:
1154
1155            rows:   list of row numbers to be flagged. Default is no row
1156                    (must be explicitly specified to execute row-based
1157                    flagging).
1158
1159            unflag: if True, unflag the data.
1160
1161        """
1162        varlist = vars()
1163        if rows is None:
1164            rows = []
1165        self._flag_row(rows, unflag)
1166        self._add_history("flag_row", varlist)
1167
1168    @asaplog_post_dec
1169    def clip(self, uthres=None, dthres=None, clipoutside=True, unflag=False):
1170        """\
1171        Flag the selected data outside a specified range (in channel-base)
1172
1173        Parameters:
1174
1175            uthres:      upper threshold.
1176
1177            dthres:      lower threshold
1178
1179            clipoutside: True for flagging data outside the range
1180                         [dthres:uthres].
1181                         False for flagging data inside the range.
1182
1183            unflag:      if True, unflag the data.
1184
1185        """
1186        varlist = vars()
1187        self._clip(uthres, dthres, clipoutside, unflag)
1188        self._add_history("clip", varlist)
1189
1190    @asaplog_post_dec
1191    def lag_flag(self, start, end, unit="MHz", insitu=None):
1192        """\
1193        Flag the data in 'lag' space by providing a frequency to remove.
1194        Flagged data in the scantable get interpolated over the region.
1195        No taper is applied.
1196
1197        Parameters:
1198
1199            start:    the start frequency (really a period within the
1200                      bandwidth)  or period to remove
1201
1202            end:      the end frequency or period to remove
1203
1204            unit:     the frequency unit (default "MHz") or "" for
1205                      explicit lag channels
1206
1207        *Notes*:
1208
1209            It is recommended to flag edges of the band or strong
1210            signals beforehand.
1211
1212        """
1213        if insitu is None: insitu = rcParams['insitu']
1214        self._math._setinsitu(insitu)
1215        varlist = vars()
1216        base = { "GHz": 1000000000., "MHz": 1000000., "kHz": 1000., "Hz": 1.}
1217        if not (unit == "" or base.has_key(unit)):
1218            raise ValueError("%s is not a valid unit." % unit)
1219        if unit == "":
1220            s = scantable(self._math._lag_flag(self, start, end, "lags"))
1221        else:
1222            s = scantable(self._math._lag_flag(self, start*base[unit],
1223                                               end*base[unit], "frequency"))
1224        s._add_history("lag_flag", varlist)
1225        if insitu:
1226            self._assign(s)
1227        else:
1228            return s
1229
1230    @asaplog_post_dec
1231    def fft(self, rowno=None, mask=None, getrealimag=False):
1232        """\
1233        Apply FFT to the spectra.
1234        Flagged data in the scantable get interpolated over the region.
1235
1236        Parameters:
1237
1238            rowno:          The row number(s) to be processed. int, list
1239                            and tuple are accepted. By default (None), FFT
1240                            is applied to the whole data.
1241
1242            mask:           Auxiliary channel mask(s). Given as a boolean
1243                            list, it is applied to all specified rows.
1244                            A list of boolean lists can also be used to
1245                            apply different masks. In the latter case, the
1246                            length of 'mask' must be the same as that of
1247                            'rowno'. The default is None.
1248       
1249            getrealimag:    If True, returns the real and imaginary part
1250                            values of the complex results.
1251                            If False (the default), returns the amplitude
1252                            (absolute value) normalised with Ndata/2 and
1253                            phase (argument, in unit of radian).
1254
1255        Returns:
1256
1257            A list of dictionaries containing the results for each spectrum.
1258            Each dictionary contains two values, the real and the imaginary
1259            parts when getrealimag = True, or the amplitude(absolute value)
1260            and the phase(argument) when getrealimag = False. The key for
1261            these values are 'real' and 'imag', or 'ampl' and 'phase',
1262            respectively.
1263        """
1264        if rowno is None:
1265            rowno = []
1266        if isinstance(rowno, int):
1267            rowno = [rowno]
1268        elif not (isinstance(rowno, list) or isinstance(rowno, tuple)):
1269            raise TypeError("The row number(s) must be int, list or tuple.")
1270        if len(rowno) == 0: rowno = [i for i in xrange(self.nrow())]
1271
1272        usecommonmask = True
1273       
1274        if mask is None:
1275            mask = []
1276        if isinstance(mask, list) or isinstance(mask, tuple):
1277            if len(mask) == 0:
1278                mask = [[]]
1279            else:
1280                if isinstance(mask[0], bool):
1281                    if len(mask) != self.nchan(self.getif(rowno[0])):
1282                        raise ValueError("The spectra and the mask have "
1283                                         "different length.")
1284                    mask = [mask]
1285                elif isinstance(mask[0], list) or isinstance(mask[0], tuple):
1286                    usecommonmask = False
1287                    if len(mask) != len(rowno):
1288                        raise ValueError("When specifying masks for each "
1289                                         "spectrum, the numbers of them "
1290                                         "must be identical.")
1291                    for i in xrange(mask):
1292                        if len(mask[i]) != self.nchan(self.getif(rowno[i])):
1293                            raise ValueError("The spectra and the mask have "
1294                                             "different length.")
1295                else:
1296                    raise TypeError("The mask must be a boolean list or "
1297                                    "a list of boolean list.")
1298        else:
1299            raise TypeError("The mask must be a boolean list or a list of "
1300                            "boolean list.")
1301
1302        res = []
1303
1304        imask = 0
1305        for whichrow in rowno:
1306            fspec = self._fft(whichrow, mask[imask], getrealimag)
1307            nspec = len(fspec)
1308           
1309            i = 0
1310            v1 = []
1311            v2 = []
1312            reselem = {"real":[],"imag":[]} if getrealimag \
1313                                            else {"ampl":[],"phase":[]}
1314           
1315            while (i < nspec):
1316                v1.append(fspec[i])
1317                v2.append(fspec[i+1])
1318                i += 2
1319           
1320            if getrealimag:
1321                reselem["real"]  += v1
1322                reselem["imag"]  += v2
1323            else:
1324                reselem["ampl"]  += v1
1325                reselem["phase"] += v2
1326           
1327            res.append(reselem)
1328           
1329            if not usecommonmask:
1330                imask += 1
1331       
1332        return res
1333
1334    @asaplog_post_dec
1335    def create_mask(self, *args, **kwargs):
1336        """\
1337        Compute and return a mask based on [min, max] windows.
1338        The specified windows are to be INCLUDED, when the mask is
1339        applied.
1340
1341        Parameters:
1342
1343            [min, max], [min2, max2], ...
1344                Pairs of start/end points (inclusive)specifying the regions
1345                to be masked
1346
1347            invert:     optional argument. If specified as True,
1348                        return an inverted mask, i.e. the regions
1349                        specified are EXCLUDED
1350
1351            row:        create the mask using the specified row for
1352                        unit conversions, default is row=0
1353                        only necessary if frequency varies over rows.
1354
1355        Examples::
1356
1357            scan.set_unit('channel')
1358            # a)
1359            msk = scan.create_mask([400, 500], [800, 900])
1360            # masks everything outside 400 and 500
1361            # and 800 and 900 in the unit 'channel'
1362
1363            # b)
1364            msk = scan.create_mask([400, 500], [800, 900], invert=True)
1365            # masks the regions between 400 and 500
1366            # and 800 and 900 in the unit 'channel'
1367
1368            # c)
1369            #mask only channel 400
1370            msk =  scan.create_mask([400])
1371
1372        """
1373        row = kwargs.get("row", 0)
1374        data = self._getabcissa(row)
1375        u = self._getcoordinfo()[0]
1376        if u == "":
1377            u = "channel"
1378        msg = "The current mask window unit is %s" % u
1379        i = self._check_ifs()
1380        if not i:
1381            msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1382        asaplog.push(msg)
1383        n = len(data)
1384        msk = _n_bools(n, False)
1385        # test if args is a 'list' or a 'normal *args - UGLY!!!
1386
1387        ws = (isinstance(args[-1][-1], int)
1388              or isinstance(args[-1][-1], float)) and args or args[0]
1389        for window in ws:
1390            if len(window) == 1:
1391                window = [window[0], window[0]]
1392            if len(window) == 0 or  len(window) > 2:
1393                raise ValueError("A window needs to be defined as "
1394                                 "[start(, end)]")
1395            if window[0] > window[1]:
1396                tmp = window[0]
1397                window[0] = window[1]
1398                window[1] = tmp
1399            for i in range(n):
1400                if data[i] >= window[0] and data[i] <= window[1]:
1401                    msk[i] = True
1402        if kwargs.has_key('invert'):
1403            if kwargs.get('invert'):
1404                msk = mask_not(msk)
1405        return msk
1406
1407    def get_masklist(self, mask=None, row=0, silent=False):
1408        """\
1409        Compute and return a list of mask windows, [min, max].
1410
1411        Parameters:
1412
1413            mask:       channel mask, created with create_mask.
1414
1415            row:        calcutate the masklist using the specified row
1416                        for unit conversions, default is row=0
1417                        only necessary if frequency varies over rows.
1418
1419        Returns:
1420
1421            [min, max], [min2, max2], ...
1422                Pairs of start/end points (inclusive)specifying
1423                the masked regions
1424
1425        """
1426        if not (isinstance(mask,list) or isinstance(mask, tuple)):
1427            raise TypeError("The mask should be list or tuple.")
1428        if len(mask) <= 0:
1429            raise TypeError("The mask elements should be > 0")
1430        data = self._getabcissa(row)
1431        if len(data) != len(mask):
1432            msg = "Number of channels in scantable != number of mask elements"
1433            raise TypeError(msg)
1434        u = self._getcoordinfo()[0]
1435        if u == "":
1436            u = "channel"
1437        msg = "The current mask window unit is %s" % u
1438        i = self._check_ifs()
1439        if not i:
1440            msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1441        if not silent:
1442            asaplog.push(msg)
1443        masklist = []
1444        ist, ien = None, None
1445        ist, ien=self.get_mask_indices(mask)
1446        if ist is not None and ien is not None:
1447            for i in xrange(len(ist)):
1448                range=[data[ist[i]],data[ien[i]]]
1449                range.sort()
1450                masklist.append([range[0],range[1]])
1451        return masklist
1452
1453    def get_mask_indices(self, mask=None):
1454        """\
1455        Compute and Return lists of mask start indices and mask end indices.
1456
1457        Parameters:
1458
1459            mask:       channel mask, created with create_mask.
1460
1461        Returns:
1462
1463            List of mask start indices and that of mask end indices,
1464            i.e., [istart1,istart2,....], [iend1,iend2,....].
1465
1466        """
1467        if not (isinstance(mask,list) or isinstance(mask, tuple)):
1468            raise TypeError("The mask should be list or tuple.")
1469        if len(mask) <= 0:
1470            raise TypeError("The mask elements should be > 0")
1471        istart = []
1472        iend = []
1473        if mask[0]:
1474            istart.append(0)
1475        for i in range(len(mask)-1):
1476            if not mask[i] and mask[i+1]:
1477                istart.append(i+1)
1478            elif mask[i] and not mask[i+1]:
1479                iend.append(i)
1480        if mask[len(mask)-1]:
1481            iend.append(len(mask)-1)
1482        if len(istart) != len(iend):
1483            raise RuntimeError("Numbers of mask start != mask end.")
1484        for i in range(len(istart)):
1485            if istart[i] > iend[i]:
1486                raise RuntimeError("Mask start index > mask end index")
1487                break
1488        return istart,iend
1489
1490    @asaplog_post_dec
1491    def parse_maskexpr(self, maskstring):
1492        """
1493        Parse CASA type mask selection syntax (IF dependent).
1494
1495        Parameters:
1496            maskstring : A string mask selection expression.
1497                         A comma separated selections mean different IF -
1498                         channel combinations. IFs and channel selections
1499                         are partitioned by a colon, ':'.
1500                     examples:
1501                         ''          = all IFs (all channels)
1502                         '<2,4~6,9'  = IFs 0,1,4,5,6,9 (all channels)
1503                         '3:3~45;60' = channels 3 to 45 and 60 in IF 3
1504                         '0~1:2~6,8' = channels 2 to 6 in IFs 0,1, and
1505                                       all channels in IF8
1506        Returns:
1507        A dictionary of selected (valid) IF and masklist pairs,
1508        e.g. {'0': [[50,250],[350,462]], '2': [[100,400],[550,974]]}
1509        """
1510        if not isinstance(maskstring,str):
1511            asaplog.post()
1512            asaplog.push("Invalid mask expression")
1513            asaplog.post("ERROR")
1514       
1515        valid_ifs = self.getifnos()
1516        frequnit = self.get_unit()
1517        seldict = {}
1518        if maskstring == "":
1519            maskstring = str(valid_ifs)[1:-1]
1520        ## split each selection
1521        sellist = maskstring.split(',')
1522        for currselstr in sellist:
1523            selset = currselstr.split(':')
1524            # spw and mask string (may include ~, < or >)
1525            spwmasklist = self._parse_selection(selset[0], typestr='integer',
1526                                                offset=1, minval=min(valid_ifs),
1527                                                maxval=max(valid_ifs))
1528            for spwlist in spwmasklist:
1529                selspws = []
1530                for ispw in range(spwlist[0],spwlist[1]+1):
1531                    # Put into the list only if ispw exists
1532                    if valid_ifs.count(ispw):
1533                        selspws.append(ispw)
1534            del spwmasklist, spwlist
1535
1536            # parse frequency mask list
1537            if len(selset) > 1:
1538                freqmasklist = self._parse_selection(selset[1], typestr='float',
1539                                                     offset=0.)
1540            else:
1541                # want to select the whole spectrum
1542                freqmasklist = [None]
1543
1544            ## define a dictionary of spw - masklist combination
1545            for ispw in selspws:
1546                #print "working on", ispw
1547                spwstr = str(ispw)
1548                if len(selspws) == 0:
1549                    # empty spw
1550                    continue
1551                else:
1552                    ## want to get min and max of the spw and
1553                    ## offset to set for '<' and '>'
1554                    if frequnit == 'channel':
1555                        minfreq = 0
1556                        maxfreq = self.nchan(ifno=ispw)
1557                        offset = 0.5
1558                    else:
1559                        ## This is ugly part. need improvement
1560                        for ifrow in xrange(self.nrow()):
1561                            if self.getif(ifrow) == ispw:
1562                                #print "IF",ispw,"found in row =",ifrow
1563                                break
1564                        freqcoord = self.get_coordinate(ifrow)
1565                        freqs = self._getabcissa(ifrow)
1566                        minfreq = min(freqs)
1567                        maxfreq = max(freqs)
1568                        if len(freqs) == 1:
1569                            offset = 0.5
1570                        elif frequnit.find('Hz') > 0:
1571                            offset = abs(freqcoord.to_frequency(1,
1572                                                                unit=frequnit)
1573                                         -freqcoord.to_frequency(0,
1574                                                                 unit=frequnit)
1575                                         )*0.5
1576                        elif frequnit.find('m/s') > 0:
1577                            offset = abs(freqcoord.to_velocity(1,
1578                                                               unit=frequnit)
1579                                         -freqcoord.to_velocity(0,
1580                                                                unit=frequnit)
1581                                         )*0.5
1582                        else:
1583                            asaplog.post()
1584                            asaplog.push("Invalid frequency unit")
1585                            asaplog.post("ERROR")
1586                        del freqs, freqcoord, ifrow
1587                    for freq in freqmasklist:
1588                        selmask = freq or [minfreq, maxfreq]
1589                        if selmask[0] == None:
1590                            ## selection was "<freq[1]".
1591                            if selmask[1] < minfreq:
1592                                ## avoid adding region selection
1593                                selmask = None
1594                            else:
1595                                selmask = [minfreq,selmask[1]-offset]
1596                        elif selmask[1] == None:
1597                            ## selection was ">freq[0]"
1598                            if selmask[0] > maxfreq:
1599                                ## avoid adding region selection
1600                                selmask = None
1601                            else:
1602                                selmask = [selmask[0]+offset,maxfreq]
1603                        if selmask:
1604                            if not seldict.has_key(spwstr):
1605                                # new spw selection
1606                                seldict[spwstr] = []
1607                            seldict[spwstr] += [selmask]
1608                    del minfreq,maxfreq,offset,freq,selmask
1609                del spwstr
1610            del freqmasklist
1611        del valid_ifs
1612        if len(seldict) == 0:
1613            asaplog.post()
1614            asaplog.push("No valid selection in the mask expression: "
1615                         +maskstring)
1616            asaplog.post("WARN")
1617            return None
1618        msg = "Selected masklist:\n"
1619        for sif, lmask in seldict.iteritems():
1620            msg += "   IF"+sif+" - "+str(lmask)+"\n"
1621        asaplog.push(msg)
1622        return seldict
1623
1624    def _parse_selection(self, selstr, typestr='float', offset=0.,
1625                         minval=None, maxval=None):
1626        """
1627        Parameters:
1628            selstr :    The Selection string, e.g., '<3;5~7;100~103;9'
1629            typestr :   The type of the values in returned list
1630                        ('integer' or 'float')
1631            offset :    The offset value to subtract from or add to
1632                        the boundary value if the selection string
1633                        includes '<' or '>'
1634            minval, maxval :  The minimum/maximum values to set if the
1635                              selection string includes '<' or '>'.
1636                              The list element is filled with None by default.
1637        Returns:
1638            A list of min/max pair of selections.
1639        Example:
1640            _parseSelection('<3;5~7;9',typestr='int',offset=1,minval=0)
1641            returns [[0,2],[5,7],[9,9]]
1642        """
1643        selgroups = selstr.split(';')
1644        sellists = []
1645        if typestr.lower().startswith('int'):
1646            formatfunc = int
1647        else:
1648            formatfunc = float
1649       
1650        for currsel in  selgroups:
1651            if currsel.find('~') > 0:
1652                minsel = formatfunc(currsel.split('~')[0].strip())
1653                maxsel = formatfunc(currsel.split('~')[1].strip())
1654            elif currsel.strip().startswith('<'):
1655                minsel = minval
1656                maxsel = formatfunc(currsel.split('<')[1].strip()) \
1657                         - formatfunc(offset)
1658            elif currsel.strip().startswith('>'):
1659                minsel = formatfunc(currsel.split('>')[1].strip()) \
1660                         + formatfunc(offset)
1661                maxsel = maxval
1662            else:
1663                minsel = formatfunc(currsel)
1664                maxsel = formatfunc(currsel)
1665            sellists.append([minsel,maxsel])
1666        return sellists
1667
1668#    def get_restfreqs(self):
1669#        """
1670#        Get the restfrequency(s) stored in this scantable.
1671#        The return value(s) are always of unit 'Hz'
1672#        Parameters:
1673#            none
1674#        Returns:
1675#            a list of doubles
1676#        """
1677#        return list(self._getrestfreqs())
1678
1679    def get_restfreqs(self, ids=None):
1680        """\
1681        Get the restfrequency(s) stored in this scantable.
1682        The return value(s) are always of unit 'Hz'
1683
1684        Parameters:
1685
1686            ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
1687                 be retrieved
1688
1689        Returns:
1690
1691            dictionary containing ids and a list of doubles for each id
1692
1693        """
1694        if ids is None:
1695            rfreqs = {}
1696            idlist = self.getmolnos()
1697            for i in idlist:
1698                rfreqs[i] = list(self._getrestfreqs(i))
1699            return rfreqs
1700        else:
1701            if type(ids) == list or type(ids) == tuple:
1702                rfreqs = {}
1703                for i in ids:
1704                    rfreqs[i] = list(self._getrestfreqs(i))
1705                return rfreqs
1706            else:
1707                return list(self._getrestfreqs(ids))
1708
1709    @asaplog_post_dec
1710    def set_restfreqs(self, freqs=None, unit='Hz'):
1711        """\
1712        Set or replace the restfrequency specified and
1713        if the 'freqs' argument holds a scalar,
1714        then that rest frequency will be applied to all the selected
1715        data.  If the 'freqs' argument holds
1716        a vector, then it MUST be of equal or smaller length than
1717        the number of IFs (and the available restfrequencies will be
1718        replaced by this vector).  In this case, *all* data have
1719        the restfrequency set per IF according
1720        to the corresponding value you give in the 'freqs' vector.
1721        E.g. 'freqs=[1e9, 2e9]'  would mean IF 0 gets restfreq 1e9 and
1722        IF 1 gets restfreq 2e9.
1723
1724        You can also specify the frequencies via a linecatalog.
1725
1726        Parameters:
1727
1728            freqs:   list of rest frequency values or string idenitfiers
1729
1730            unit:    unit for rest frequency (default 'Hz')
1731
1732
1733        Example::
1734
1735            # set the given restfrequency for the all currently selected IFs
1736            scan.set_restfreqs(freqs=1.4e9)
1737            # set restfrequencies for the n IFs  (n > 1) in the order of the
1738            # list, i.e
1739            # IF0 -> 1.4e9, IF1 ->  1.41e9, IF3 -> 1.42e9
1740            # len(list_of_restfreqs) == nIF
1741            # for nIF == 1 the following will set multiple restfrequency for
1742            # that IF
1743            scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1744            # set multiple restfrequencies per IF. as a list of lists where
1745            # the outer list has nIF elements, the inner s arbitrary
1746            scan.set_restfreqs(freqs=[[1.4e9, 1.41e9], [1.67e9]])
1747
1748       *Note*:
1749
1750            To do more sophisticate Restfrequency setting, e.g. on a
1751            source and IF basis, use scantable.set_selection() before using
1752            this function::
1753
1754                # provided your scantable is called scan
1755                selection = selector()
1756                selection.set_name("ORION*")
1757                selection.set_ifs([1])
1758                scan.set_selection(selection)
1759                scan.set_restfreqs(freqs=86.6e9)
1760
1761        """
1762        varlist = vars()
1763        from asap import linecatalog
1764        # simple  value
1765        if isinstance(freqs, int) or isinstance(freqs, float):
1766            self._setrestfreqs([freqs], [""], unit)
1767        # list of values
1768        elif isinstance(freqs, list) or isinstance(freqs, tuple):
1769            # list values are scalars
1770            if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1771                if len(freqs) == 1:
1772                    self._setrestfreqs(freqs, [""], unit)
1773                else:
1774                    # allow the 'old' mode of setting mulitple IFs
1775                    sel = selector()
1776                    savesel = self._getselection()
1777                    iflist = self.getifnos()
1778                    if len(freqs)>len(iflist):
1779                        raise ValueError("number of elements in list of list "
1780                                         "exeeds the current IF selections")
1781                    iflist = self.getifnos()
1782                    for i, fval in enumerate(freqs):
1783                        sel.set_ifs(iflist[i])
1784                        self._setselection(sel)
1785                        self._setrestfreqs([fval], [""], unit)
1786                    self._setselection(savesel)
1787
1788            # list values are dict, {'value'=, 'name'=)
1789            elif isinstance(freqs[-1], dict):
1790                values = []
1791                names = []
1792                for d in freqs:
1793                    values.append(d["value"])
1794                    names.append(d["name"])
1795                self._setrestfreqs(values, names, unit)
1796            elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1797                sel = selector()
1798                savesel = self._getselection()
1799                iflist = self.getifnos()
1800                if len(freqs)>len(iflist):
1801                    raise ValueError("number of elements in list of list exeeds"
1802                                     " the current IF selections")
1803                for i, fval in enumerate(freqs):
1804                    sel.set_ifs(iflist[i])
1805                    self._setselection(sel)
1806                    self._setrestfreqs(fval, [""], unit)
1807                self._setselection(savesel)
1808        # freqs are to be taken from a linecatalog
1809        elif isinstance(freqs, linecatalog):
1810            sel = selector()
1811            savesel = self._getselection()
1812            for i in xrange(freqs.nrow()):
1813                sel.set_ifs(iflist[i])
1814                self._setselection(sel)
1815                self._setrestfreqs([freqs.get_frequency(i)],
1816                                   [freqs.get_name(i)], "MHz")
1817                # ensure that we are not iterating past nIF
1818                if i == self.nif()-1: break
1819            self._setselection(savesel)
1820        else:
1821            return
1822        self._add_history("set_restfreqs", varlist)
1823
1824    @asaplog_post_dec
1825    def shift_refpix(self, delta):
1826        """\
1827        Shift the reference pixel of the Spectra Coordinate by an
1828        integer amount.
1829
1830        Parameters:
1831
1832            delta:   the amount to shift by
1833
1834        *Note*:
1835
1836            Be careful using this with broadband data.
1837
1838        """
1839        varlist = vars()
1840        Scantable.shift_refpix(self, delta)
1841        s._add_history("shift_refpix", varlist)
1842
1843    @asaplog_post_dec
1844    def history(self, filename=None):
1845        """\
1846        Print the history. Optionally to a file.
1847
1848        Parameters:
1849
1850            filename:    The name of the file to save the history to.
1851
1852        """
1853        hist = list(self._gethistory())
1854        out = "-"*80
1855        for h in hist:
1856            if h.startswith("---"):
1857                out = "\n".join([out, h])
1858            else:
1859                items = h.split("##")
1860                date = items[0]
1861                func = items[1]
1862                items = items[2:]
1863                out += "\n"+date+"\n"
1864                out += "Function: %s\n  Parameters:" % (func)
1865                for i in items:
1866                    if i == '':
1867                        continue
1868                    s = i.split("=")
1869                    out += "\n   %s = %s" % (s[0], s[1])
1870                out = "\n".join([out, "-"*80])
1871        if filename is not None:
1872            if filename is "":
1873                filename = 'scantable_history.txt'
1874            import os
1875            filename = os.path.expandvars(os.path.expanduser(filename))
1876            if not os.path.isdir(filename):
1877                data = open(filename, 'w')
1878                data.write(out)
1879                data.close()
1880            else:
1881                msg = "Illegal file name '%s'." % (filename)
1882                raise IOError(msg)
1883        return page(out)
1884
1885    #
1886    # Maths business
1887    #
1888    @asaplog_post_dec
1889    def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1890        """\
1891        Return the (time) weighted average of a scan. Scans will be averaged
1892        only if the source direction (RA/DEC) is within 1' otherwise
1893
1894        *Note*:
1895
1896            in channels only - align if necessary
1897
1898        Parameters:
1899
1900            mask:     an optional mask (only used for 'var' and 'tsys'
1901                      weighting)
1902
1903            scanav:   True averages each scan separately
1904                      False (default) averages all scans together,
1905
1906            weight:   Weighting scheme.
1907                      'none'     (mean no weight)
1908                      'var'      (1/var(spec) weighted)
1909                      'tsys'     (1/Tsys**2 weighted)
1910                      'tint'     (integration time weighted)
1911                      'tintsys'  (Tint/Tsys**2)
1912                      'median'   ( median averaging)
1913                      The default is 'tint'
1914
1915            align:    align the spectra in velocity before averaging. It takes
1916                      the time of the first spectrum as reference time.
1917
1918        Example::
1919
1920            # time average the scantable without using a mask
1921            newscan = scan.average_time()
1922
1923        """
1924        varlist = vars()
1925        weight = weight or 'TINT'
1926        mask = mask or ()
1927        scanav = (scanav and 'SCAN') or 'NONE'
1928        scan = (self, )
1929
1930        if align:
1931            scan = (self.freq_align(insitu=False), )
1932        s = None
1933        if weight.upper() == 'MEDIAN':
1934            s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1935                                                     scanav))
1936        else:
1937            s = scantable(self._math._average(scan, mask, weight.upper(),
1938                          scanav))
1939        s._add_history("average_time", varlist)
1940        return s
1941
1942    @asaplog_post_dec
1943    def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
1944        """\
1945        Return a scan where all spectra are converted to either
1946        Jansky or Kelvin depending upon the flux units of the scan table.
1947        By default the function tries to look the values up internally.
1948        If it can't find them (or if you want to over-ride), you must
1949        specify EITHER jyperk OR eta (and D which it will try to look up
1950        also if you don't set it). jyperk takes precedence if you set both.
1951
1952        Parameters:
1953
1954            jyperk:      the Jy / K conversion factor
1955
1956            eta:         the aperture efficiency
1957
1958            d:           the geometric diameter (metres)
1959
1960            insitu:      if False a new scantable is returned.
1961                         Otherwise, the scaling is done in-situ
1962                         The default is taken from .asaprc (False)
1963
1964        """
1965        if insitu is None: insitu = rcParams['insitu']
1966        self._math._setinsitu(insitu)
1967        varlist = vars()
1968        jyperk = jyperk or -1.0
1969        d = d or -1.0
1970        eta = eta or -1.0
1971        s = scantable(self._math._convertflux(self, d, eta, jyperk))
1972        s._add_history("convert_flux", varlist)
1973        if insitu: self._assign(s)
1974        else: return s
1975
1976    @asaplog_post_dec
1977    def gain_el(self, poly=None, filename="", method="linear", insitu=None):
1978        """\
1979        Return a scan after applying a gain-elevation correction.
1980        The correction can be made via either a polynomial or a
1981        table-based interpolation (and extrapolation if necessary).
1982        You specify polynomial coefficients, an ascii table or neither.
1983        If you specify neither, then a polynomial correction will be made
1984        with built in coefficients known for certain telescopes (an error
1985        will occur if the instrument is not known).
1986        The data and Tsys are *divided* by the scaling factors.
1987
1988        Parameters:
1989
1990            poly:        Polynomial coefficients (default None) to compute a
1991                         gain-elevation correction as a function of
1992                         elevation (in degrees).
1993
1994            filename:    The name of an ascii file holding correction factors.
1995                         The first row of the ascii file must give the column
1996                         names and these MUST include columns
1997                         "ELEVATION" (degrees) and "FACTOR" (multiply data
1998                         by this) somewhere.
1999                         The second row must give the data type of the
2000                         column. Use 'R' for Real and 'I' for Integer.
2001                         An example file would be
2002                         (actual factors are arbitrary) :
2003
2004                         TIME ELEVATION FACTOR
2005                         R R R
2006                         0.1 0 0.8
2007                         0.2 20 0.85
2008                         0.3 40 0.9
2009                         0.4 60 0.85
2010                         0.5 80 0.8
2011                         0.6 90 0.75
2012
2013            method:      Interpolation method when correcting from a table.
2014                         Values are  "nearest", "linear" (default), "cubic"
2015                         and "spline"
2016
2017            insitu:      if False a new scantable is returned.
2018                         Otherwise, the scaling is done in-situ
2019                         The default is taken from .asaprc (False)
2020
2021        """
2022
2023        if insitu is None: insitu = rcParams['insitu']
2024        self._math._setinsitu(insitu)
2025        varlist = vars()
2026        poly = poly or ()
2027        from os.path import expandvars
2028        filename = expandvars(filename)
2029        s = scantable(self._math._gainel(self, poly, filename, method))
2030        s._add_history("gain_el", varlist)
2031        if insitu:
2032            self._assign(s)
2033        else:
2034            return s
2035
2036    @asaplog_post_dec
2037    def freq_align(self, reftime=None, method='cubic', insitu=None):
2038        """\
2039        Return a scan where all rows have been aligned in frequency/velocity.
2040        The alignment frequency frame (e.g. LSRK) is that set by function
2041        set_freqframe.
2042
2043        Parameters:
2044
2045            reftime:     reference time to align at. By default, the time of
2046                         the first row of data is used.
2047
2048            method:      Interpolation method for regridding the spectra.
2049                         Choose from "nearest", "linear", "cubic" (default)
2050                         and "spline"
2051
2052            insitu:      if False a new scantable is returned.
2053                         Otherwise, the scaling is done in-situ
2054                         The default is taken from .asaprc (False)
2055
2056        """
2057        if insitu is None: insitu = rcParams["insitu"]
2058        oldInsitu = self._math._insitu()
2059        self._math._setinsitu(insitu)
2060        varlist = vars()
2061        reftime = reftime or ""
2062        s = scantable(self._math._freq_align(self, reftime, method))
2063        s._add_history("freq_align", varlist)
2064        self._math._setinsitu(oldInsitu)
2065        if insitu:
2066            self._assign(s)
2067        else:
2068            return s
2069
2070    @asaplog_post_dec
2071    def opacity(self, tau=None, insitu=None):
2072        """\
2073        Apply an opacity correction. The data
2074        and Tsys are multiplied by the correction factor.
2075
2076        Parameters:
2077
2078            tau:         (list of) opacity from which the correction factor is
2079                         exp(tau*ZD)
2080                         where ZD is the zenith-distance.
2081                         If a list is provided, it has to be of length nIF,
2082                         nIF*nPol or 1 and in order of IF/POL, e.g.
2083                         [opif0pol0, opif0pol1, opif1pol0 ...]
2084                         if tau is `None` the opacities are determined from a
2085                         model.
2086
2087            insitu:      if False a new scantable is returned.
2088                         Otherwise, the scaling is done in-situ
2089                         The default is taken from .asaprc (False)
2090
2091        """
2092        if insitu is None:
2093            insitu = rcParams['insitu']
2094        self._math._setinsitu(insitu)
2095        varlist = vars()
2096        if not hasattr(tau, "__len__"):
2097            tau = [tau]
2098        s = scantable(self._math._opacity(self, tau))
2099        s._add_history("opacity", varlist)
2100        if insitu:
2101            self._assign(s)
2102        else:
2103            return s
2104
2105    @asaplog_post_dec
2106    def bin(self, width=5, insitu=None):
2107        """\
2108        Return a scan where all spectra have been binned up.
2109
2110        Parameters:
2111
2112            width:       The bin width (default=5) in pixels
2113
2114            insitu:      if False a new scantable is returned.
2115                         Otherwise, the scaling is done in-situ
2116                         The default is taken from .asaprc (False)
2117
2118        """
2119        if insitu is None:
2120            insitu = rcParams['insitu']
2121        self._math._setinsitu(insitu)
2122        varlist = vars()
2123        s = scantable(self._math._bin(self, width))
2124        s._add_history("bin", varlist)
2125        if insitu:
2126            self._assign(s)
2127        else:
2128            return s
2129
2130    @asaplog_post_dec
2131    def resample(self, width=5, method='cubic', insitu=None):
2132        """\
2133        Return a scan where all spectra have been binned up.
2134
2135        Parameters:
2136
2137            width:       The bin width (default=5) in pixels
2138
2139            method:      Interpolation method when correcting from a table.
2140                         Values are  "nearest", "linear", "cubic" (default)
2141                         and "spline"
2142
2143            insitu:      if False a new scantable is returned.
2144                         Otherwise, the scaling is done in-situ
2145                         The default is taken from .asaprc (False)
2146
2147        """
2148        if insitu is None:
2149            insitu = rcParams['insitu']
2150        self._math._setinsitu(insitu)
2151        varlist = vars()
2152        s = scantable(self._math._resample(self, method, width))
2153        s._add_history("resample", varlist)
2154        if insitu:
2155            self._assign(s)
2156        else:
2157            return s
2158
2159    @asaplog_post_dec
2160    def average_pol(self, mask=None, weight='none'):
2161        """\
2162        Average the Polarisations together.
2163
2164        Parameters:
2165
2166            mask:        An optional mask defining the region, where the
2167                         averaging will be applied. The output will have all
2168                         specified points masked.
2169
2170            weight:      Weighting scheme. 'none' (default), 'var' (1/var(spec)
2171                         weighted), or 'tsys' (1/Tsys**2 weighted)
2172
2173        """
2174        varlist = vars()
2175        mask = mask or ()
2176        s = scantable(self._math._averagepol(self, mask, weight.upper()))
2177        s._add_history("average_pol", varlist)
2178        return s
2179
2180    @asaplog_post_dec
2181    def average_beam(self, mask=None, weight='none'):
2182        """\
2183        Average the Beams together.
2184
2185        Parameters:
2186            mask:        An optional mask defining the region, where the
2187                         averaging will be applied. The output will have all
2188                         specified points masked.
2189
2190            weight:      Weighting scheme. 'none' (default), 'var' (1/var(spec)
2191                         weighted), or 'tsys' (1/Tsys**2 weighted)
2192
2193        """
2194        varlist = vars()
2195        mask = mask or ()
2196        s = scantable(self._math._averagebeams(self, mask, weight.upper()))
2197        s._add_history("average_beam", varlist)
2198        return s
2199
2200    def parallactify(self, pflag):
2201        """\
2202        Set a flag to indicate whether this data should be treated as having
2203        been 'parallactified' (total phase == 0.0)
2204
2205        Parameters:
2206
2207            pflag:  Bool indicating whether to turn this on (True) or
2208                    off (False)
2209
2210        """
2211        varlist = vars()
2212        self._parallactify(pflag)
2213        self._add_history("parallactify", varlist)
2214
2215    @asaplog_post_dec
2216    def convert_pol(self, poltype=None):
2217        """\
2218        Convert the data to a different polarisation type.
2219        Note that you will need cross-polarisation terms for most conversions.
2220
2221        Parameters:
2222
2223            poltype:    The new polarisation type. Valid types are:
2224                        "linear", "circular", "stokes" and "linpol"
2225
2226        """
2227        varlist = vars()
2228        s = scantable(self._math._convertpol(self, poltype))
2229        s._add_history("convert_pol", varlist)
2230        return s
2231
2232    @asaplog_post_dec
2233    def smooth(self, kernel="hanning", width=5.0, order=2, plot=False,
2234               insitu=None):
2235        """\
2236        Smooth the spectrum by the specified kernel (conserving flux).
2237
2238        Parameters:
2239
2240            kernel:     The type of smoothing kernel. Select from
2241                        'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
2242                        or 'poly'
2243
2244            width:      The width of the kernel in pixels. For hanning this is
2245                        ignored otherwise it defauls to 5 pixels.
2246                        For 'gaussian' it is the Full Width Half
2247                        Maximum. For 'boxcar' it is the full width.
2248                        For 'rmedian' and 'poly' it is the half width.
2249
2250            order:      Optional parameter for 'poly' kernel (default is 2), to
2251                        specify the order of the polnomial. Ignored by all other
2252                        kernels.
2253
2254            plot:       plot the original and the smoothed spectra.
2255                        In this each indivual fit has to be approved, by
2256                        typing 'y' or 'n'
2257
2258            insitu:     if False a new scantable is returned.
2259                        Otherwise, the scaling is done in-situ
2260                        The default is taken from .asaprc (False)
2261
2262        """
2263        if insitu is None: insitu = rcParams['insitu']
2264        self._math._setinsitu(insitu)
2265        varlist = vars()
2266
2267        if plot: orgscan = self.copy()
2268
2269        s = scantable(self._math._smooth(self, kernel.lower(), width, order))
2270        s._add_history("smooth", varlist)
2271
2272        if plot:
2273            from asap.asapplotter import new_asaplot
2274            theplot = new_asaplot(rcParams['plotter.gui'])
2275            theplot.set_panels()
2276            ylab=s._get_ordinate_label()
2277            #theplot.palette(0,["#777777","red"])
2278            for r in xrange(s.nrow()):
2279                xsm=s._getabcissa(r)
2280                ysm=s._getspectrum(r)
2281                xorg=orgscan._getabcissa(r)
2282                yorg=orgscan._getspectrum(r)
2283                theplot.clear()
2284                theplot.hold()
2285                theplot.set_axes('ylabel',ylab)
2286                theplot.set_axes('xlabel',s._getabcissalabel(r))
2287                theplot.set_axes('title',s._getsourcename(r))
2288                theplot.set_line(label='Original',color="#777777")
2289                theplot.plot(xorg,yorg)
2290                theplot.set_line(label='Smoothed',color="red")
2291                theplot.plot(xsm,ysm)
2292                ### Ugly part for legend
2293                for i in [0,1]:
2294                    theplot.subplots[0]['lines'].append(
2295                        [theplot.subplots[0]['axes'].lines[i]]
2296                        )
2297                theplot.release()
2298                ### Ugly part for legend
2299                theplot.subplots[0]['lines']=[]
2300                res = raw_input("Accept smoothing ([y]/n): ")
2301                if res.upper() == 'N':
2302                    s._setspectrum(yorg, r)
2303            theplot.quit()
2304            del theplot
2305            del orgscan
2306
2307        if insitu: self._assign(s)
2308        else: return s
2309
2310    @asaplog_post_dec
2311    def _parse_wn(self, wn):
2312        if isinstance(wn, list) or isinstance(wn, tuple):
2313            return wn
2314        elif isinstance(wn, int):
2315            return [ wn ]
2316        elif isinstance(wn, str):
2317            if '-' in wn:                            # case 'a-b' : return [a,a+1,...,b-1,b]
2318                val = wn.split('-')
2319                val = [int(val[0]), int(val[1])]
2320                val.sort()
2321                res = [i for i in xrange(val[0], val[1]+1)]
2322            elif wn[:2] == '<=' or wn[:2] == '=<':   # cases '<=a','=<a' : return [0,1,...,a-1,a]
2323                val = int(wn[2:])+1
2324                res = [i for i in xrange(val)]
2325            elif wn[-2:] == '>=' or wn[-2:] == '=>': # cases 'a>=','a=>' : return [0,1,...,a-1,a]
2326                val = int(wn[:-2])+1
2327                res = [i for i in xrange(val)]
2328            elif wn[0] == '<':                       # case '<a' :         return [0,1,...,a-2,a-1]
2329                val = int(wn[1:])
2330                res = [i for i in xrange(val)]
2331            elif wn[-1] == '>':                      # case 'a>' :         return [0,1,...,a-2,a-1]
2332                val = int(wn[:-1])
2333                res = [i for i in xrange(val)]
2334            elif wn[:2] == '>=' or wn[:2] == '=>':   # cases '>=a','=>a' : return [a,-999], which is
2335                                                     #                     then interpreted in C++
2336                                                     #                     side as [a,a+1,...,a_nyq]
2337                                                     #                     (CAS-3759)
2338                val = int(wn[2:])
2339                res = [val, -999]
2340                #res = [i for i in xrange(val, self.nchan()/2+1)]
2341            elif wn[-2:] == '<=' or wn[-2:] == '=<': # cases 'a<=','a=<' : return [a,-999], which is
2342                                                     #                     then interpreted in C++
2343                                                     #                     side as [a,a+1,...,a_nyq]
2344                                                     #                     (CAS-3759)
2345                val = int(wn[:-2])
2346                res = [val, -999]
2347                #res = [i for i in xrange(val, self.nchan()/2+1)]
2348            elif wn[0] == '>':                       # case '>a' :         return [a+1,-999], which is
2349                                                     #                     then interpreted in C++
2350                                                     #                     side as [a+1,a+2,...,a_nyq]
2351                                                     #                     (CAS-3759)
2352                val = int(wn[1:])+1
2353                res = [val, -999]
2354                #res = [i for i in xrange(val, self.nchan()/2+1)]
2355            elif wn[-1] == '<':                      # case 'a<' :         return [a+1,-999], which is
2356                                                     #                     then interpreted in C++
2357                                                     #                     side as [a+1,a+2,...,a_nyq]
2358                                                     #                     (CAS-3759)
2359                val = int(wn[:-1])+1
2360                res = [val, -999]
2361                #res = [i for i in xrange(val, self.nchan()/2+1)]
2362
2363            return res
2364        else:
2365            msg = 'wrong value given for addwn/rejwn'
2366            raise RuntimeError(msg)
2367
2368
2369    @asaplog_post_dec
2370    def sinusoid_baseline(self, insitu=None, mask=None, applyfft=None,
2371                          fftmethod=None, fftthresh=None,
2372                          addwn=None, rejwn=None, clipthresh=None,
2373                          clipniter=None, plot=None,
2374                          getresidual=None, showprogress=None,
2375                          minnrow=None, outlog=None, blfile=None):
2376        """\
2377        Return a scan which has been baselined (all rows) with sinusoidal
2378        functions.
2379
2380        Parameters:
2381            insitu:        if False a new scantable is returned.
2382                           Otherwise, the scaling is done in-situ
2383                           The default is taken from .asaprc (False)
2384            mask:          an optional mask
2385            applyfft:      if True use some method, such as FFT, to find
2386                           strongest sinusoidal components in the wavenumber
2387                           domain to be used for baseline fitting.
2388                           default is True.
2389            fftmethod:     method to find the strong sinusoidal components.
2390                           now only 'fft' is available and it is the default.
2391            fftthresh:     the threshold to select wave numbers to be used for
2392                           fitting from the distribution of amplitudes in the
2393                           wavenumber domain.
2394                           both float and string values accepted.
2395                           given a float value, the unit is set to sigma.
2396                           for string values, allowed formats include:
2397                               'xsigma' or 'x' (= x-sigma level. e.g.,
2398                               '3sigma'), or
2399                               'topx' (= the x strongest ones, e.g. 'top5').
2400                           default is 3.0 (unit: sigma).
2401            addwn:         the additional wave numbers to be used for fitting.
2402                           list or integer value is accepted to specify every
2403                           wave numbers. also string value can be used in case
2404                           you need to specify wave numbers in a certain range,
2405                           e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
2406                                 '<a'  (= 0,1,...,a-2,a-1),
2407                                 '>=a' (= a, a+1, ... up to the maximum wave
2408                                        number corresponding to the Nyquist
2409                                        frequency for the case of FFT).
2410                           default is [0].
2411            rejwn:         the wave numbers NOT to be used for fitting.
2412                           can be set just as addwn but has higher priority:
2413                           wave numbers which are specified both in addwn
2414                           and rejwn will NOT be used. default is [].
2415            clipthresh:    Clipping threshold. (default is 3.0, unit: sigma)
2416            clipniter:     maximum number of iteration of 'clipthresh'-sigma
2417                           clipping (default is 0)
2418            plot:      *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2419                           plot the fit and the residual. In this each
2420                           indivual fit has to be approved, by typing 'y'
2421                           or 'n'
2422            getresidual:   if False, returns best-fit values instead of
2423                           residual. (default is True)
2424            showprogress:  show progress status for large data.
2425                           default is True.
2426            minnrow:       minimum number of input spectra to show.
2427                           default is 1000.
2428            outlog:        Output the coefficients of the best-fit
2429                           function to logger (default is False)
2430            blfile:        Name of a text file in which the best-fit
2431                           parameter values to be written
2432                           (default is '': no file/logger output)
2433
2434        Example:
2435            # return a scan baselined by a combination of sinusoidal curves
2436            # having wave numbers in spectral window up to 10,
2437            # also with 3-sigma clipping, iteration up to 4 times
2438            bscan = scan.sinusoid_baseline(addwn='<=10',clipthresh=3.0,clipniter=4)
2439
2440        Note:
2441            The best-fit parameter values output in logger and/or blfile are now
2442            based on specunit of 'channel'.
2443        """
2444       
2445        try:
2446            varlist = vars()
2447       
2448            if insitu is None: insitu = rcParams['insitu']
2449            if insitu:
2450                workscan = self
2451            else:
2452                workscan = self.copy()
2453           
2454            #if mask          is None: mask          = [True for i in xrange(workscan.nchan())]
2455            if mask          is None: mask          = []
2456            if applyfft      is None: applyfft      = True
2457            if fftmethod     is None: fftmethod     = 'fft'
2458            if fftthresh     is None: fftthresh     = 3.0
2459            if addwn         is None: addwn         = [0]
2460            if rejwn         is None: rejwn         = []
2461            if clipthresh    is None: clipthresh    = 3.0
2462            if clipniter     is None: clipniter     = 0
2463            if plot          is None: plot          = False
2464            if getresidual   is None: getresidual   = True
2465            if showprogress  is None: showprogress  = True
2466            if minnrow       is None: minnrow       = 1000
2467            if outlog        is None: outlog        = False
2468            if blfile        is None: blfile        = ''
2469
2470            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2471            workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(),
2472                                        str(fftthresh).lower(),
2473                                        workscan._parse_wn(addwn),
2474                                        workscan._parse_wn(rejwn), clipthresh,
2475                                        clipniter, getresidual,
2476                                        pack_progress_params(showprogress,
2477                                                             minnrow), outlog,
2478                                        blfile)
2479            workscan._add_history('sinusoid_baseline', varlist)
2480           
2481            if insitu:
2482                self._assign(workscan)
2483            else:
2484                return workscan
2485           
2486        except RuntimeError, e:
2487            raise_fitting_failure_exception(e)
2488
2489
2490    @asaplog_post_dec
2491    def auto_sinusoid_baseline(self, insitu=None, mask=None, applyfft=None,
2492                               fftmethod=None, fftthresh=None,
2493                               addwn=None, rejwn=None, clipthresh=None,
2494                               clipniter=None, edge=None, threshold=None,
2495                               chan_avg_limit=None, plot=None,
2496                               getresidual=None, showprogress=None,
2497                               minnrow=None, outlog=None, blfile=None):
2498        """\
2499        Return a scan which has been baselined (all rows) with sinusoidal
2500        functions.
2501        Spectral lines are detected first using linefinder and masked out
2502        to avoid them affecting the baseline solution.
2503
2504        Parameters:
2505            insitu:         if False a new scantable is returned.
2506                            Otherwise, the scaling is done in-situ
2507                            The default is taken from .asaprc (False)
2508            mask:           an optional mask retreived from scantable
2509            applyfft:       if True use some method, such as FFT, to find
2510                            strongest sinusoidal components in the wavenumber
2511                            domain to be used for baseline fitting.
2512                            default is True.
2513            fftmethod:      method to find the strong sinusoidal components.
2514                            now only 'fft' is available and it is the default.
2515            fftthresh:      the threshold to select wave numbers to be used for
2516                            fitting from the distribution of amplitudes in the
2517                            wavenumber domain.
2518                            both float and string values accepted.
2519                            given a float value, the unit is set to sigma.
2520                            for string values, allowed formats include:
2521                                'xsigma' or 'x' (= x-sigma level. e.g.,
2522                                '3sigma'), or
2523                                'topx' (= the x strongest ones, e.g. 'top5').
2524                            default is 3.0 (unit: sigma).
2525            addwn:          the additional wave numbers to be used for fitting.
2526                            list or integer value is accepted to specify every
2527                            wave numbers. also string value can be used in case
2528                            you need to specify wave numbers in a certain range,
2529                            e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
2530                                  '<a'  (= 0,1,...,a-2,a-1),
2531                                  '>=a' (= a, a+1, ... up to the maximum wave
2532                                         number corresponding to the Nyquist
2533                                         frequency for the case of FFT).
2534                            default is [0].
2535            rejwn:          the wave numbers NOT to be used for fitting.
2536                            can be set just as addwn but has higher priority:
2537                            wave numbers which are specified both in addwn
2538                            and rejwn will NOT be used. default is [].
2539            clipthresh:     Clipping threshold. (default is 3.0, unit: sigma)
2540            clipniter:      maximum number of iteration of 'clipthresh'-sigma
2541                            clipping (default is 0)
2542            edge:           an optional number of channel to drop at
2543                            the edge of spectrum. If only one value is
2544                            specified, the same number will be dropped
2545                            from both sides of the spectrum. Default
2546                            is to keep all channels. Nested tuples
2547                            represent individual edge selection for
2548                            different IFs (a number of spectral channels
2549                            can be different)
2550            threshold:      the threshold used by line finder. It is
2551                            better to keep it large as only strong lines
2552                            affect the baseline solution.
2553            chan_avg_limit: a maximum number of consequtive spectral
2554                            channels to average during the search of
2555                            weak and broad lines. The default is no
2556                            averaging (and no search for weak lines).
2557                            If such lines can affect the fitted baseline
2558                            (e.g. a high order polynomial is fitted),
2559                            increase this parameter (usually values up
2560                            to 8 are reasonable). Most users of this
2561                            method should find the default value sufficient.
2562            plot:       *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2563                            plot the fit and the residual. In this each
2564                            indivual fit has to be approved, by typing 'y'
2565                            or 'n'
2566            getresidual:    if False, returns best-fit values instead of
2567                            residual. (default is True)
2568            showprogress:   show progress status for large data.
2569                            default is True.
2570            minnrow:        minimum number of input spectra to show.
2571                            default is 1000.
2572            outlog:         Output the coefficients of the best-fit
2573                            function to logger (default is False)
2574            blfile:         Name of a text file in which the best-fit
2575                            parameter values to be written
2576                            (default is "": no file/logger output)
2577
2578        Example:
2579            bscan = scan.auto_sinusoid_baseline(addwn='<=10', insitu=False)
2580       
2581        Note:
2582            The best-fit parameter values output in logger and/or blfile are now
2583            based on specunit of 'channel'.
2584        """
2585
2586        try:
2587            varlist = vars()
2588
2589            if insitu is None: insitu = rcParams['insitu']
2590            if insitu:
2591                workscan = self
2592            else:
2593                workscan = self.copy()
2594           
2595            #if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
2596            if mask           is None: mask           = []
2597            if applyfft       is None: applyfft       = True
2598            if fftmethod      is None: fftmethod      = 'fft'
2599            if fftthresh      is None: fftthresh      = 3.0
2600            if addwn          is None: addwn          = [0]
2601            if rejwn          is None: rejwn          = []
2602            if clipthresh     is None: clipthresh     = 3.0
2603            if clipniter      is None: clipniter      = 0
2604            if edge           is None: edge           = (0,0)
2605            if threshold      is None: threshold      = 3
2606            if chan_avg_limit is None: chan_avg_limit = 1
2607            if plot           is None: plot           = False
2608            if getresidual    is None: getresidual    = True
2609            if showprogress   is None: showprogress   = True
2610            if minnrow        is None: minnrow        = 1000
2611            if outlog         is None: outlog         = False
2612            if blfile         is None: blfile         = ''
2613
2614            #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2615            workscan._auto_sinusoid_baseline(mask, applyfft,
2616                                             fftmethod.lower(),
2617                                             str(fftthresh).lower(),
2618                                             workscan._parse_wn(addwn),
2619                                             workscan._parse_wn(rejwn),
2620                                             clipthresh, clipniter,
2621                                             normalise_edge_param(edge),
2622                                             threshold, chan_avg_limit,
2623                                             getresidual,
2624                                             pack_progress_params(showprogress,
2625                                                                  minnrow),
2626                                             outlog, blfile)
2627            workscan._add_history("auto_sinusoid_baseline", varlist)
2628           
2629            if insitu:
2630                self._assign(workscan)
2631            else:
2632                return workscan
2633           
2634        except RuntimeError, e:
2635            raise_fitting_failure_exception(e)
2636
2637    @asaplog_post_dec
2638    def cspline_baseline(self, insitu=None, mask=None, npiece=None,
2639                         clipthresh=None, clipniter=None, plot=None,
2640                         getresidual=None, showprogress=None, minnrow=None,
2641                         outlog=None, blfile=None):
2642        """\
2643        Return a scan which has been baselined (all rows) by cubic spline
2644        function (piecewise cubic polynomial).
2645
2646        Parameters:
2647            insitu:       If False a new scantable is returned.
2648                          Otherwise, the scaling is done in-situ
2649                          The default is taken from .asaprc (False)
2650            mask:         An optional mask
2651            npiece:       Number of pieces. (default is 2)
2652            clipthresh:   Clipping threshold. (default is 3.0, unit: sigma)
2653            clipniter:    maximum number of iteration of 'clipthresh'-sigma
2654                          clipping (default is 0)
2655            plot:     *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2656                          plot the fit and the residual. In this each
2657                          indivual fit has to be approved, by typing 'y'
2658                          or 'n'
2659            getresidual:  if False, returns best-fit values instead of
2660                          residual. (default is True)
2661            showprogress: show progress status for large data.
2662                          default is True.
2663            minnrow:      minimum number of input spectra to show.
2664                          default is 1000.
2665            outlog:       Output the coefficients of the best-fit
2666                          function to logger (default is False)
2667            blfile:       Name of a text file in which the best-fit
2668                          parameter values to be written
2669                          (default is "": no file/logger output)
2670
2671        Example:
2672            # return a scan baselined by a cubic spline consisting of 2 pieces
2673            # (i.e., 1 internal knot),
2674            # also with 3-sigma clipping, iteration up to 4 times
2675            bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
2676       
2677        Note:
2678            The best-fit parameter values output in logger and/or blfile are now
2679            based on specunit of 'channel'.
2680        """
2681       
2682        try:
2683            varlist = vars()
2684           
2685            if insitu is None: insitu = rcParams['insitu']
2686            if insitu:
2687                workscan = self
2688            else:
2689                workscan = self.copy()
2690
2691            #if mask         is None: mask         = [True for i in xrange(workscan.nchan())]
2692            if mask         is None: mask         = []
2693            if npiece       is None: npiece       = 2
2694            if clipthresh   is None: clipthresh   = 3.0
2695            if clipniter    is None: clipniter    = 0
2696            if plot         is None: plot         = False
2697            if getresidual  is None: getresidual  = True
2698            if showprogress is None: showprogress = True
2699            if minnrow      is None: minnrow      = 1000
2700            if outlog       is None: outlog       = False
2701            if blfile       is None: blfile       = ''
2702
2703            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2704            workscan._cspline_baseline(mask, npiece, clipthresh, clipniter,
2705                                       getresidual,
2706                                       pack_progress_params(showprogress,
2707                                                            minnrow), outlog,
2708                                       blfile)
2709            workscan._add_history("cspline_baseline", varlist)
2710           
2711            if insitu:
2712                self._assign(workscan)
2713            else:
2714                return workscan
2715           
2716        except RuntimeError, e:
2717            raise_fitting_failure_exception(e)
2718
2719    @asaplog_post_dec
2720    def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None,
2721                              clipthresh=None, clipniter=None,
2722                              edge=None, threshold=None, chan_avg_limit=None,
2723                              getresidual=None, plot=None,
2724                              showprogress=None, minnrow=None, outlog=None,
2725                              blfile=None):
2726        """\
2727        Return a scan which has been baselined (all rows) by cubic spline
2728        function (piecewise cubic polynomial).
2729        Spectral lines are detected first using linefinder and masked out
2730        to avoid them affecting the baseline solution.
2731
2732        Parameters:
2733            insitu:         if False a new scantable is returned.
2734                            Otherwise, the scaling is done in-situ
2735                            The default is taken from .asaprc (False)
2736            mask:           an optional mask retreived from scantable
2737            npiece:         Number of pieces. (default is 2)
2738            clipthresh:     Clipping threshold. (default is 3.0, unit: sigma)
2739            clipniter:      maximum number of iteration of 'clipthresh'-sigma
2740                            clipping (default is 0)
2741            edge:           an optional number of channel to drop at
2742                            the edge of spectrum. If only one value is
2743                            specified, the same number will be dropped
2744                            from both sides of the spectrum. Default
2745                            is to keep all channels. Nested tuples
2746                            represent individual edge selection for
2747                            different IFs (a number of spectral channels
2748                            can be different)
2749            threshold:      the threshold used by line finder. It is
2750                            better to keep it large as only strong lines
2751                            affect the baseline solution.
2752            chan_avg_limit: a maximum number of consequtive spectral
2753                            channels to average during the search of
2754                            weak and broad lines. The default is no
2755                            averaging (and no search for weak lines).
2756                            If such lines can affect the fitted baseline
2757                            (e.g. a high order polynomial is fitted),
2758                            increase this parameter (usually values up
2759                            to 8 are reasonable). Most users of this
2760                            method should find the default value sufficient.
2761            plot:       *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2762                            plot the fit and the residual. In this each
2763                            indivual fit has to be approved, by typing 'y'
2764                            or 'n'
2765            getresidual:    if False, returns best-fit values instead of
2766                            residual. (default is True)
2767            showprogress:   show progress status for large data.
2768                            default is True.
2769            minnrow:        minimum number of input spectra to show.
2770                            default is 1000.
2771            outlog:         Output the coefficients of the best-fit
2772                            function to logger (default is False)
2773            blfile:         Name of a text file in which the best-fit
2774                            parameter values to be written
2775                            (default is "": no file/logger output)
2776
2777        Example:
2778            bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
2779       
2780        Note:
2781            The best-fit parameter values output in logger and/or blfile are now
2782            based on specunit of 'channel'.
2783        """
2784
2785        try:
2786            varlist = vars()
2787
2788            if insitu is None: insitu = rcParams['insitu']
2789            if insitu:
2790                workscan = self
2791            else:
2792                workscan = self.copy()
2793           
2794            #if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
2795            if mask           is None: mask           = []
2796            if npiece         is None: npiece         = 2
2797            if clipthresh     is None: clipthresh     = 3.0
2798            if clipniter      is None: clipniter      = 0
2799            if edge           is None: edge           = (0, 0)
2800            if threshold      is None: threshold      = 3
2801            if chan_avg_limit is None: chan_avg_limit = 1
2802            if plot           is None: plot           = False
2803            if getresidual    is None: getresidual    = True
2804            if showprogress   is None: showprogress   = True
2805            if minnrow        is None: minnrow        = 1000
2806            if outlog         is None: outlog         = False
2807            if blfile         is None: blfile         = ''
2808
2809            #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2810            workscan._auto_cspline_baseline(mask, npiece, clipthresh,
2811                                            clipniter,
2812                                            normalise_edge_param(edge),
2813                                            threshold,
2814                                            chan_avg_limit, getresidual,
2815                                            pack_progress_params(showprogress,
2816                                                                 minnrow),
2817                                            outlog, blfile)
2818            workscan._add_history("auto_cspline_baseline", varlist)
2819           
2820            if insitu:
2821                self._assign(workscan)
2822            else:
2823                return workscan
2824           
2825        except RuntimeError, e:
2826            raise_fitting_failure_exception(e)
2827
2828    @asaplog_post_dec
2829    def poly_baseline(self, mask=None, order=None, insitu=None, plot=None,
2830                      getresidual=None, showprogress=None, minnrow=None,
2831                      outlog=None, blfile=None):
2832        """\
2833        Return a scan which has been baselined (all rows) by a polynomial.
2834        Parameters:
2835            insitu:       if False a new scantable is returned.
2836                          Otherwise, the scaling is done in-situ
2837                          The default is taken from .asaprc (False)
2838            mask:         an optional mask
2839            order:        the order of the polynomial (default is 0)
2840            plot:         plot the fit and the residual. In this each
2841                          indivual fit has to be approved, by typing 'y'
2842                          or 'n'
2843            getresidual:  if False, returns best-fit values instead of
2844                          residual. (default is True)
2845            showprogress: show progress status for large data.
2846                          default is True.
2847            minnrow:      minimum number of input spectra to show.
2848                          default is 1000.
2849            outlog:       Output the coefficients of the best-fit
2850                          function to logger (default is False)
2851            blfile:       Name of a text file in which the best-fit
2852                          parameter values to be written
2853                          (default is "": no file/logger output)
2854
2855        Example:
2856            # return a scan baselined by a third order polynomial,
2857            # not using a mask
2858            bscan = scan.poly_baseline(order=3)
2859        """
2860       
2861        try:
2862            varlist = vars()
2863       
2864            if insitu is None:
2865                insitu = rcParams["insitu"]
2866            if insitu:
2867                workscan = self
2868            else:
2869                workscan = self.copy()
2870
2871            #if mask         is None: mask         = [True for i in \
2872            #                                           xrange(workscan.nchan())]
2873            if mask         is None: mask         = []
2874            if order        is None: order        = 0
2875            if plot         is None: plot         = False
2876            if getresidual  is None: getresidual  = True
2877            if showprogress is None: showprogress = True
2878            if minnrow      is None: minnrow      = 1000
2879            if outlog       is None: outlog       = False
2880            if blfile       is None: blfile       = ""
2881
2882            if plot:
2883                outblfile = (blfile != "") and \
2884                    os.path.exists(os.path.expanduser(
2885                                    os.path.expandvars(blfile))
2886                                   )
2887                if outblfile:
2888                    blf = open(blfile, "a")
2889               
2890                f = fitter()
2891                f.set_function(lpoly=order)
2892               
2893                rows = xrange(workscan.nrow())
2894                #if len(rows) > 0: workscan._init_blinfo()
2895               
2896                for r in rows:
2897                    f.x = workscan._getabcissa(r)
2898                    f.y = workscan._getspectrum(r)
2899                    f.mask = mask_and(mask, workscan._getmask(r))    # (CAS-1434)
2900                    f.data = None
2901                    f.fit()
2902                   
2903                    f.plot(residual=True)
2904                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
2905                    if accept_fit.upper() == "N":
2906                        #workscan._append_blinfo(None, None, None)
2907                        continue
2908                   
2909                    blpars = f.get_parameters()
2910                    masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2911                    #workscan._append_blinfo(blpars, masklist, f.mask)
2912                    workscan._setspectrum((f.fitter.getresidual()
2913                                           if getresidual else f.fitter.getfit()), r)
2914                   
2915                    if outblfile:
2916                        rms = workscan.get_rms(f.mask, r)
2917                        dataout = \
2918                            workscan.format_blparams_row(blpars["params"],
2919                                                         blpars["fixed"],
2920                                                         rms, str(masklist),
2921                                                         r, True)
2922                        blf.write(dataout)
2923
2924                f._p.unmap()
2925                f._p = None
2926
2927                if outblfile:
2928                    blf.close()
2929            else:
2930                workscan._poly_baseline(mask, order, getresidual,
2931                                        pack_progress_params(showprogress,
2932                                                             minnrow),
2933                                        outlog, blfile)
2934           
2935            workscan._add_history("poly_baseline", varlist)
2936           
2937            if insitu:
2938                self._assign(workscan)
2939            else:
2940                return workscan
2941           
2942        except RuntimeError, e:
2943            raise_fitting_failure_exception(e)
2944
2945    @asaplog_post_dec
2946    def auto_poly_baseline(self, mask=None, order=None, edge=None,
2947                           threshold=None, chan_avg_limit=None,
2948                           plot=None, insitu=None,
2949                           getresidual=None, showprogress=None,
2950                           minnrow=None, outlog=None, blfile=None):
2951        """\
2952        Return a scan which has been baselined (all rows) by a polynomial.
2953        Spectral lines are detected first using linefinder and masked out
2954        to avoid them affecting the baseline solution.
2955
2956        Parameters:
2957            insitu:         if False a new scantable is returned.
2958                            Otherwise, the scaling is done in-situ
2959                            The default is taken from .asaprc (False)
2960            mask:           an optional mask retreived from scantable
2961            order:          the order of the polynomial (default is 0)
2962            edge:           an optional number of channel to drop at
2963                            the edge of spectrum. If only one value is
2964                            specified, the same number will be dropped
2965                            from both sides of the spectrum. Default
2966                            is to keep all channels. Nested tuples
2967                            represent individual edge selection for
2968                            different IFs (a number of spectral channels
2969                            can be different)
2970            threshold:      the threshold used by line finder. It is
2971                            better to keep it large as only strong lines
2972                            affect the baseline solution.
2973            chan_avg_limit: a maximum number of consequtive spectral
2974                            channels to average during the search of
2975                            weak and broad lines. The default is no
2976                            averaging (and no search for weak lines).
2977                            If such lines can affect the fitted baseline
2978                            (e.g. a high order polynomial is fitted),
2979                            increase this parameter (usually values up
2980                            to 8 are reasonable). Most users of this
2981                            method should find the default value sufficient.
2982            plot:           plot the fit and the residual. In this each
2983                            indivual fit has to be approved, by typing 'y'
2984                            or 'n'
2985            getresidual:    if False, returns best-fit values instead of
2986                            residual. (default is True)
2987            showprogress:   show progress status for large data.
2988                            default is True.
2989            minnrow:        minimum number of input spectra to show.
2990                            default is 1000.
2991            outlog:         Output the coefficients of the best-fit
2992                            function to logger (default is False)
2993            blfile:         Name of a text file in which the best-fit
2994                            parameter values to be written
2995                            (default is "": no file/logger output)
2996
2997        Example:
2998            bscan = scan.auto_poly_baseline(order=7, insitu=False)
2999        """
3000
3001        try:
3002            varlist = vars()
3003
3004            if insitu is None:
3005                insitu = rcParams['insitu']
3006            if insitu:
3007                workscan = self
3008            else:
3009                workscan = self.copy()
3010
3011            #if mask           is None: mask           = [True for i in xrange(workscan.nchan())]
3012            if mask           is None: mask           = []
3013            if order          is None: order          = 0
3014            if edge           is None: edge           = (0, 0)
3015            if threshold      is None: threshold      = 3
3016            if chan_avg_limit is None: chan_avg_limit = 1
3017            if plot           is None: plot           = False
3018            if getresidual    is None: getresidual    = True
3019            if showprogress   is None: showprogress   = True
3020            if minnrow        is None: minnrow        = 1000
3021            if outlog         is None: outlog         = False
3022            if blfile         is None: blfile         = ''
3023
3024            edge = normalise_edge_param(edge)
3025
3026            if plot:
3027                outblfile = (blfile != "") and \
3028                    os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
3029                if outblfile: blf = open(blfile, "a")
3030               
3031                from asap.asaplinefind import linefinder
3032                fl = linefinder()
3033                fl.set_options(threshold=threshold, avg_limit=chan_avg_limit)
3034                fl.set_scan(workscan)
3035               
3036                f = fitter()
3037                f.set_function(lpoly=order)
3038
3039                rows = xrange(workscan.nrow())
3040                #if len(rows) > 0: workscan._init_blinfo()
3041               
3042                for r in rows:
3043                    idx = 2*workscan.getif(r)
3044                    fl.find_lines(r, mask_and(mask, workscan._getmask(r)),
3045                                  edge[idx:idx+2])  # (CAS-1434)
3046
3047                    f.x = workscan._getabcissa(r)
3048                    f.y = workscan._getspectrum(r)
3049                    f.mask = fl.get_mask()
3050                    f.data = None
3051                    f.fit()
3052
3053                    f.plot(residual=True)
3054                    accept_fit = raw_input("Accept fit ( [y]/n ): ")
3055                    if accept_fit.upper() == "N":
3056                        #workscan._append_blinfo(None, None, None)
3057                        continue
3058
3059                    blpars = f.get_parameters()
3060                    masklist = workscan.get_masklist(f.mask, row=r, silent=True)
3061                    #workscan._append_blinfo(blpars, masklist, f.mask)
3062                    workscan._setspectrum(
3063                        (f.fitter.getresidual() if getresidual
3064                                                else f.fitter.getfit()), r
3065                        )
3066
3067                    if outblfile:
3068                        rms = workscan.get_rms(f.mask, r)
3069                        dataout = \
3070                            workscan.format_blparams_row(blpars["params"],
3071                                                         blpars["fixed"],
3072                                                         rms, str(masklist),
3073                                                         r, True)
3074                        blf.write(dataout)
3075                   
3076                f._p.unmap()
3077                f._p = None
3078
3079                if outblfile: blf.close()
3080            else:
3081                workscan._auto_poly_baseline(mask, order, edge, threshold,
3082                                             chan_avg_limit, getresidual,
3083                                             pack_progress_params(showprogress,
3084                                                                  minnrow),
3085                                             outlog, blfile)
3086
3087            workscan._add_history("auto_poly_baseline", varlist)
3088           
3089            if insitu:
3090                self._assign(workscan)
3091            else:
3092                return workscan
3093           
3094        except RuntimeError, e:
3095            raise_fitting_failure_exception(e)
3096
3097    def _init_blinfo(self):
3098        """\
3099        Initialise the following three auxiliary members:
3100           blpars : parameters of the best-fit baseline,
3101           masklists : mask data (edge positions of masked channels) and
3102           actualmask : mask data (in boolean list),
3103        to keep for use later (including output to logger/text files).
3104        Used by poly_baseline() and auto_poly_baseline() in case of
3105        'plot=True'.
3106        """
3107        self.blpars = []
3108        self.masklists = []
3109        self.actualmask = []
3110        return
3111
3112    def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
3113        """\
3114        Append baseline-fitting related info to blpars, masklist and
3115        actualmask.
3116        """
3117        self.blpars.append(data_blpars)
3118        self.masklists.append(data_masklists)
3119        self.actualmask.append(data_actualmask)
3120        return
3121       
3122    @asaplog_post_dec
3123    def rotate_linpolphase(self, angle):
3124        """\
3125        Rotate the phase of the complex polarization O=Q+iU correlation.
3126        This is always done in situ in the raw data.  So if you call this
3127        function more than once then each call rotates the phase further.
3128
3129        Parameters:
3130
3131            angle:   The angle (degrees) to rotate (add) by.
3132
3133        Example::
3134
3135            scan.rotate_linpolphase(2.3)
3136
3137        """
3138        varlist = vars()
3139        self._math._rotate_linpolphase(self, angle)
3140        self._add_history("rotate_linpolphase", varlist)
3141        return
3142
3143    @asaplog_post_dec
3144    def rotate_xyphase(self, angle):
3145        """\
3146        Rotate the phase of the XY correlation.  This is always done in situ
3147        in the data.  So if you call this function more than once
3148        then each call rotates the phase further.
3149
3150        Parameters:
3151
3152            angle:   The angle (degrees) to rotate (add) by.
3153
3154        Example::
3155
3156            scan.rotate_xyphase(2.3)
3157
3158        """
3159        varlist = vars()
3160        self._math._rotate_xyphase(self, angle)
3161        self._add_history("rotate_xyphase", varlist)
3162        return
3163
3164    @asaplog_post_dec
3165    def swap_linears(self):
3166        """\
3167        Swap the linear polarisations XX and YY, or better the first two
3168        polarisations as this also works for ciculars.
3169        """
3170        varlist = vars()
3171        self._math._swap_linears(self)
3172        self._add_history("swap_linears", varlist)
3173        return
3174
3175    @asaplog_post_dec
3176    def invert_phase(self):
3177        """\
3178        Invert the phase of the complex polarisation
3179        """
3180        varlist = vars()
3181        self._math._invert_phase(self)
3182        self._add_history("invert_phase", varlist)
3183        return
3184
3185    @asaplog_post_dec
3186    def add(self, offset, insitu=None):
3187        """\
3188        Return a scan where all spectra have the offset added
3189
3190        Parameters:
3191
3192            offset:      the offset
3193
3194            insitu:      if False a new scantable is returned.
3195                         Otherwise, the scaling is done in-situ
3196                         The default is taken from .asaprc (False)
3197
3198        """
3199        if insitu is None: insitu = rcParams['insitu']
3200        self._math._setinsitu(insitu)
3201        varlist = vars()
3202        s = scantable(self._math._unaryop(self, offset, "ADD", False))
3203        s._add_history("add", varlist)
3204        if insitu:
3205            self._assign(s)
3206        else:
3207            return s
3208
3209    @asaplog_post_dec
3210    def scale(self, factor, tsys=True, insitu=None):
3211        """\
3212
3213        Return a scan where all spectra are scaled by the given 'factor'
3214
3215        Parameters:
3216
3217            factor:      the scaling factor (float or 1D float list)
3218
3219            insitu:      if False a new scantable is returned.
3220                         Otherwise, the scaling is done in-situ
3221                         The default is taken from .asaprc (False)
3222
3223            tsys:        if True (default) then apply the operation to Tsys
3224                         as well as the data
3225
3226        """
3227        if insitu is None: insitu = rcParams['insitu']
3228        self._math._setinsitu(insitu)
3229        varlist = vars()
3230        s = None
3231        import numpy
3232        if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
3233            if isinstance(factor[0], list) or isinstance(factor[0],
3234                                                         numpy.ndarray):
3235                from asapmath import _array2dOp
3236                s = _array2dOp( self, factor, "MUL", tsys, insitu )
3237            else:
3238                s = scantable( self._math._arrayop( self, factor,
3239                                                    "MUL", tsys ) )
3240        else:
3241            s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
3242        s._add_history("scale", varlist)
3243        if insitu:
3244            self._assign(s)
3245        else:
3246            return s
3247
3248    @preserve_selection
3249    def set_sourcetype(self, match, matchtype="pattern",
3250                       sourcetype="reference"):
3251        """\
3252        Set the type of the source to be an source or reference scan
3253        using the provided pattern.
3254
3255        Parameters:
3256
3257            match:          a Unix style pattern, regular expression or selector
3258
3259            matchtype:      'pattern' (default) UNIX style pattern or
3260                            'regex' regular expression
3261
3262            sourcetype:     the type of the source to use (source/reference)
3263
3264        """
3265        varlist = vars()
3266        basesel = self.get_selection()
3267        stype = -1
3268        if sourcetype.lower().startswith("r"):
3269            stype = 1
3270        elif sourcetype.lower().startswith("s"):
3271            stype = 0
3272        else:
3273            raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
3274        if matchtype.lower().startswith("p"):
3275            matchtype = "pattern"
3276        elif matchtype.lower().startswith("r"):
3277            matchtype = "regex"
3278        else:
3279            raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
3280        sel = selector()
3281        if isinstance(match, selector):
3282            sel = match
3283        else:
3284            sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
3285        self._setsourcetype(stype)
3286        self._add_history("set_sourcetype", varlist)
3287
3288    @asaplog_post_dec
3289    @preserve_selection
3290    def auto_quotient(self, preserve=True, mode='paired', verify=False):
3291        """\
3292        This function allows to build quotients automatically.
3293        It assumes the observation to have the same number of
3294        "ons" and "offs"
3295
3296        Parameters:
3297
3298            preserve:       you can preserve (default) the continuum or
3299                            remove it.  The equations used are
3300
3301                            preserve: Output = Toff * (on/off) - Toff
3302
3303                            remove:   Output = Toff * (on/off) - Ton
3304
3305            mode:           the on/off detection mode
3306                            'paired' (default)
3307                            identifies 'off' scans by the
3308                            trailing '_R' (Mopra/Parkes) or
3309                            '_e'/'_w' (Tid) and matches
3310                            on/off pairs from the observing pattern
3311                            'time'
3312                            finds the closest off in time
3313
3314        .. todo:: verify argument is not implemented
3315
3316        """
3317        varlist = vars()
3318        modes = ["time", "paired"]
3319        if not mode in modes:
3320            msg = "please provide valid mode. Valid modes are %s" % (modes)
3321            raise ValueError(msg)
3322        s = None
3323        if mode.lower() == "paired":
3324            sel = self.get_selection()
3325            sel.set_query("SRCTYPE==psoff")
3326            self.set_selection(sel)
3327            offs = self.copy()
3328            sel.set_query("SRCTYPE==pson")
3329            self.set_selection(sel)
3330            ons = self.copy()
3331            s = scantable(self._math._quotient(ons, offs, preserve))
3332        elif mode.lower() == "time":
3333            s = scantable(self._math._auto_quotient(self, mode, preserve))
3334        s._add_history("auto_quotient", varlist)
3335        return s
3336
3337    @asaplog_post_dec
3338    def mx_quotient(self, mask = None, weight='median', preserve=True):
3339        """\
3340        Form a quotient using "off" beams when observing in "MX" mode.
3341
3342        Parameters:
3343
3344            mask:           an optional mask to be used when weight == 'stddev'
3345
3346            weight:         How to average the off beams.  Default is 'median'.
3347
3348            preserve:       you can preserve (default) the continuum or
3349                            remove it.  The equations used are:
3350
3351                                preserve: Output = Toff * (on/off) - Toff
3352
3353                                remove:   Output = Toff * (on/off) - Ton
3354
3355        """
3356        mask = mask or ()
3357        varlist = vars()
3358        on = scantable(self._math._mx_extract(self, 'on'))
3359        preoff = scantable(self._math._mx_extract(self, 'off'))
3360        off = preoff.average_time(mask=mask, weight=weight, scanav=False)
3361        from asapmath  import quotient
3362        q = quotient(on, off, preserve)
3363        q._add_history("mx_quotient", varlist)
3364        return q
3365
3366    @asaplog_post_dec
3367    def freq_switch(self, insitu=None):
3368        """\
3369        Apply frequency switching to the data.
3370
3371        Parameters:
3372
3373            insitu:      if False a new scantable is returned.
3374                         Otherwise, the swictching is done in-situ
3375                         The default is taken from .asaprc (False)
3376
3377        """
3378        if insitu is None: insitu = rcParams['insitu']
3379        self._math._setinsitu(insitu)
3380        varlist = vars()
3381        s = scantable(self._math._freqswitch(self))
3382        s._add_history("freq_switch", varlist)
3383        if insitu:
3384            self._assign(s)
3385        else:
3386            return s
3387
3388    @asaplog_post_dec
3389    def recalc_azel(self):
3390        """Recalculate the azimuth and elevation for each position."""
3391        varlist = vars()
3392        self._recalcazel()
3393        self._add_history("recalc_azel", varlist)
3394        return
3395
3396    @asaplog_post_dec
3397    def __add__(self, other):
3398        varlist = vars()
3399        s = None
3400        if isinstance(other, scantable):
3401            s = scantable(self._math._binaryop(self, other, "ADD"))
3402        elif isinstance(other, float):
3403            s = scantable(self._math._unaryop(self, other, "ADD", False))
3404        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3405            if isinstance(other[0], list) \
3406                    or isinstance(other[0], numpy.ndarray):
3407                from asapmath import _array2dOp
3408                s = _array2dOp( self.copy(), other, "ADD", False )
3409            else:
3410                s = scantable( self._math._arrayop( self.copy(), other,
3411                                                    "ADD", False ) )
3412        else:
3413            raise TypeError("Other input is not a scantable or float value")
3414        s._add_history("operator +", varlist)
3415        return s
3416
3417    @asaplog_post_dec
3418    def __sub__(self, other):
3419        """
3420        implicit on all axes and on Tsys
3421        """
3422        varlist = vars()
3423        s = None
3424        if isinstance(other, scantable):
3425            s = scantable(self._math._binaryop(self, other, "SUB"))
3426        elif isinstance(other, float):
3427            s = scantable(self._math._unaryop(self, other, "SUB", False))
3428        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3429            if isinstance(other[0], list) \
3430                    or isinstance(other[0], numpy.ndarray):
3431                from asapmath import _array2dOp
3432                s = _array2dOp( self.copy(), other, "SUB", False )
3433            else:
3434                s = scantable( self._math._arrayop( self.copy(), other,
3435                                                    "SUB", False ) )
3436        else:
3437            raise TypeError("Other input is not a scantable or float value")
3438        s._add_history("operator -", varlist)
3439        return s
3440
3441    @asaplog_post_dec
3442    def __mul__(self, other):
3443        """
3444        implicit on all axes and on Tsys
3445        """
3446        varlist = vars()
3447        s = None
3448        if isinstance(other, scantable):
3449            s = scantable(self._math._binaryop(self, other, "MUL"))
3450        elif isinstance(other, float):
3451            s = scantable(self._math._unaryop(self, other, "MUL", False))
3452        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3453            if isinstance(other[0], list) \
3454                    or isinstance(other[0], numpy.ndarray):
3455                from asapmath import _array2dOp
3456                s = _array2dOp( self.copy(), other, "MUL", False )
3457            else:
3458                s = scantable( self._math._arrayop( self.copy(), other,
3459                                                    "MUL", False ) )
3460        else:
3461            raise TypeError("Other input is not a scantable or float value")
3462        s._add_history("operator *", varlist)
3463        return s
3464
3465
3466    @asaplog_post_dec
3467    def __div__(self, other):
3468        """
3469        implicit on all axes and on Tsys
3470        """
3471        varlist = vars()
3472        s = None
3473        if isinstance(other, scantable):
3474            s = scantable(self._math._binaryop(self, other, "DIV"))
3475        elif isinstance(other, float):
3476            if other == 0.0:
3477                raise ZeroDivisionError("Dividing by zero is not recommended")
3478            s = scantable(self._math._unaryop(self, other, "DIV", False))
3479        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3480            if isinstance(other[0], list) \
3481                    or isinstance(other[0], numpy.ndarray):
3482                from asapmath import _array2dOp
3483                s = _array2dOp( self.copy(), other, "DIV", False )
3484            else:
3485                s = scantable( self._math._arrayop( self.copy(), other,
3486                                                    "DIV", False ) )
3487        else:
3488            raise TypeError("Other input is not a scantable or float value")
3489        s._add_history("operator /", varlist)
3490        return s
3491
3492    @asaplog_post_dec
3493    def get_fit(self, row=0):
3494        """\
3495        Print or return the stored fits for a row in the scantable
3496
3497        Parameters:
3498
3499            row:    the row which the fit has been applied to.
3500
3501        """
3502        if row > self.nrow():
3503            return
3504        from asap.asapfit import asapfit
3505        fit = asapfit(self._getfit(row))
3506        asaplog.push( '%s' %(fit) )
3507        return fit.as_dict()
3508
3509    @preserve_selection
3510    def flag_nans(self):
3511        """\
3512        Utility function to flag NaN values in the scantable.
3513        """
3514        import numpy
3515        basesel = self.get_selection()
3516        for i in range(self.nrow()):
3517            sel = self.get_row_selector(i)
3518            self.set_selection(basesel+sel)
3519            nans = numpy.isnan(self._getspectrum(0))
3520        if numpy.any(nans):
3521            bnans = [ bool(v) for v in nans]
3522            self.flag(bnans)
3523
3524    def get_row_selector(self, rowno):
3525        return selector(rows=[rowno])
3526
3527    def _add_history(self, funcname, parameters):
3528        if not rcParams['scantable.history']:
3529            return
3530        # create date
3531        sep = "##"
3532        from datetime import datetime
3533        dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
3534        hist = dstr+sep
3535        hist += funcname+sep#cdate+sep
3536        if parameters.has_key('self'):
3537            del parameters['self']
3538        for k, v in parameters.iteritems():
3539            if type(v) is dict:
3540                for k2, v2 in v.iteritems():
3541                    hist += k2
3542                    hist += "="
3543                    if isinstance(v2, scantable):
3544                        hist += 'scantable'
3545                    elif k2 == 'mask':
3546                        if isinstance(v2, list) or isinstance(v2, tuple):
3547                            hist += str(self._zip_mask(v2))
3548                        else:
3549                            hist += str(v2)
3550                    else:
3551                        hist += str(v2)
3552            else:
3553                hist += k
3554                hist += "="
3555                if isinstance(v, scantable):
3556                    hist += 'scantable'
3557                elif k == 'mask':
3558                    if isinstance(v, list) or isinstance(v, tuple):
3559                        hist += str(self._zip_mask(v))
3560                    else:
3561                        hist += str(v)
3562                else:
3563                    hist += str(v)
3564            hist += sep
3565        hist = hist[:-2] # remove trailing '##'
3566        self._addhistory(hist)
3567
3568
3569    def _zip_mask(self, mask):
3570        mask = list(mask)
3571        i = 0
3572        segments = []
3573        while mask[i:].count(1):
3574            i += mask[i:].index(1)
3575            if mask[i:].count(0):
3576                j = i + mask[i:].index(0)
3577            else:
3578                j = len(mask)
3579            segments.append([i, j])
3580            i = j
3581        return segments
3582
3583    def _get_ordinate_label(self):
3584        fu = "("+self.get_fluxunit()+")"
3585        import re
3586        lbl = "Intensity"
3587        if re.match(".K.", fu):
3588            lbl = "Brightness Temperature "+ fu
3589        elif re.match(".Jy.", fu):
3590            lbl = "Flux density "+ fu
3591        return lbl
3592
3593    def _check_ifs(self):
3594#        return len(set([self.nchan(i) for i in self.getifnos()])) == 1
3595        nchans = [self.nchan(i) for i in self.getifnos()]
3596        nchans = filter(lambda t: t > 0, nchans)
3597        return (sum(nchans)/len(nchans) == nchans[0])
3598
3599    @asaplog_post_dec
3600    def _fill(self, names, unit, average, opts={}):
3601        first = True
3602        fullnames = []
3603        for name in names:
3604            name = os.path.expandvars(name)
3605            name = os.path.expanduser(name)
3606            if not os.path.exists(name):
3607                msg = "File '%s' does not exists" % (name)
3608                raise IOError(msg)
3609            fullnames.append(name)
3610        if average:
3611            asaplog.push('Auto averaging integrations')
3612        stype = int(rcParams['scantable.storage'].lower() == 'disk')
3613        for name in fullnames:
3614            tbl = Scantable(stype)
3615            if is_ms( name ):
3616                r = msfiller( tbl )
3617            else:
3618                r = filler( tbl )
3619                rx = rcParams['scantable.reference']
3620                r.setreferenceexpr(rx)
3621            msg = "Importing %s..." % (name)
3622            asaplog.push(msg, False)
3623            r.open(name, opts)
3624            r.fill()
3625            if average:
3626                tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
3627            if not first:
3628                tbl = self._math._merge([self, tbl])
3629            Scantable.__init__(self, tbl)
3630            r.close()
3631            del r, tbl
3632            first = False
3633            #flush log
3634        asaplog.post()
3635        if unit is not None:
3636            self.set_fluxunit(unit)
3637        if not is_casapy():
3638            self.set_freqframe(rcParams['scantable.freqframe'])
3639
3640
3641    def __getitem__(self, key):
3642        if key < 0:
3643            key += self.nrow()
3644        if key >= self.nrow():
3645            raise IndexError("Row index out of range.")
3646        return self._getspectrum(key)
3647
3648    def __setitem__(self, key, value):
3649        if key < 0:
3650            key += self.nrow()
3651        if key >= self.nrow():
3652            raise IndexError("Row index out of range.")
3653        if not hasattr(value, "__len__") or \
3654                len(value) > self.nchan(self.getif(key)):
3655            raise ValueError("Spectrum length doesn't match.")
3656        return self._setspectrum(value, key)
3657
3658    def __len__(self):
3659        return self.nrow()
3660
3661    def __iter__(self):
3662        for i in range(len(self)):
3663            yield self[i]
Note: See TracBrowser for help on using the repository browser.