source: trunk/python/scantable.py @ 2687

Last change on this file since 2687 was 2687, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: from asap.scantable import is_scantable

is_scantable('<any CASA image>')
# if it return False the bug is fixed

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Fixed a bug in is_scantable() that the function recognizes
CASA image as Scantable.


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