source: trunk/python/scantable.py @ 1443

Last change on this file since 1443 was 1443, checked in by Malte Marquarding, 16 years ago

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