source: trunk/python/scantable.py @ 2012

Last change on this file since 2012 was 2012, checked in by WataruKawasaki, 13 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



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