source: trunk/python/scantable.py @ 1322

Last change on this file since 1322 was 1322, checked in by mar637, 17 years ago

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