source: trunk/python/scantable.py @ 2672

Last change on this file since 2672 was 2672, checked in by Malte Marquarding, 12 years ago

introduce 'reshape' method; this includes a fix to c++ for the upper boundary assertion

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