source: trunk/python/scantable.py @ 1402

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

implemented iterator to iterate through rows

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