source: trunk/python/scantable.py @ 936

Last change on this file since 936 was 936, checked in by mar637, 18 years ago

fixed self.stm to self._math
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.6 KB
Line 
1from asap._asap import Scantable
2from asap import rcParams
3from asap import print_log
4from numarray import ones,zeros
5import sys
6
7class scantable(Scantable):
8    """
9        The ASAP container for scans
10    """
11
12    def __init__(self, filename, average=None, unit=None):
13        """
14        Create a scantable from a saved one or make a reference
15        Parameters:
16            filename:    the name of an asap table on disk
17                         or
18                         the name of a rpfits/sdfits/ms file
19                         (integrations within scans are auto averaged
20                         and the whole file is read)
21                         or
22                         [advanced] a reference to an existing
23                         scantable
24            average:     average all integrations withinb a scan on read.
25                         The default (True) is taken from .asaprc.
26            unit:         brightness unit; must be consistent with K or Jy.
27                         Over-rides the default selected by the reader
28                         (input rpfits/sdfits/ms) or replaces the value
29                         in existing scantables
30        """
31        if average is None or type(average) is not bool:
32            average = rcParams['scantable.autoaverage']
33
34        varlist = vars()
35        from asap._asap import stmath
36        self._math = stmath()
37        from asap import asaplog
38        if isinstance(filename, Scantable):
39            Scantable.__init__(self, filename)
40            if unit is not None:
41                self.set_fluxunit(unit)
42        else:
43            import os.path
44            if not os.path.exists(filename):
45                s = "File '%s' not found." % (filename)
46                if rcParams['verbose']:
47                    asaplog.push(s)
48                    print asaplog.pop().strip()
49                    return
50                raise IOError(s)
51            filename = os.path.expandvars(filename)
52            if os.path.isdir(filename):
53                # crude check if asap table
54                if os.path.exists(filename+'/table.info'):
55                    Scantable.__init__(self, filename, "memory")
56                    if unit is not None:
57                        self.set_fluxunit(unit)
58                else:
59                    msg = "The given file '%s'is not a valid asap table." % (filename)
60                    if rcParams['verbose']:
61                        print msg
62                        return
63                    else:
64                        raise IOError(msg)
65            else:
66                from asap._asap import stfiller
67                ifSel = -1
68                beamSel = -1
69                r = stfiller()
70                r._open(filename,ifSel,beamSel)
71                asaplog.push('Importing data...')
72                print_log()
73                r._read()
74                tbl = r._getdata()
75                if unit is not None:
76                    tbl.set_fluxunit(unit)
77                if average:
78                    asaplog.push('Auto averaging integrations...')
79                    print_log()
80                    tbl2 = self._math._average((tbl,),(),'NONE','SCAN', False)
81                    Scantable.__init__(self,tbl2)
82                    del tbl2
83                else:
84                    Scantable.__init__(self,tbl)
85                del r,tbl
86                self._add_history("scantable", varlist)
87        print_log()
88
89    def save(self, name=None, format=None, overwrite=False):
90        """
91        Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
92        Image FITS or MS2 format.
93        Parameters:
94            name:        the name of the outputfile. For format="FITS" this
95                         is the directory file name into which all the files
96                         will be written (default is 'asap_FITS'). For format
97                         "ASCII" this is the root file name (data in 'name'.txt
98                         and header in 'name'_header.txt)
99            format:      an optional file format. Default is ASAP.
100                         Allowed are - 'ASAP' (save as ASAP [aips++] Table),
101                                       'SDFITS' (save as SDFITS file)
102                                       'FITS' (saves each row as a FITS Image)
103                                       'ASCII' (saves as ascii text file)
104                                       'MS2' (saves as an aips++
105                                              MeasurementSet V2)
106            overwrite:   If the file should be overwritten if it exists.
107                         The default False is to return with warning
108                         without writing the output. USE WITH CARE.
109        Example:
110            scan.save('myscan.asap')
111            scan.save('myscan.sdfits','SDFITS')
112        """
113        from os import path
114        if format is None: format = rcParams['scantable.save']
115        suffix = '.'+format.lower()
116        if name is None or name =="":
117            name = 'scantable'+suffix
118            from asap import asaplog
119            msg = "No filename given. Using default name %s..." % name
120            asaplog.push(msg)
121        name = path.expandvars(name)
122        if path.isfile(name) or path.isdir(name):
123            if not overwrite:
124                msg = "File %s exists." % name
125                if rcParams['verbose']:
126                    print msg
127                    return
128                else:
129                    raise IOError(msg)
130        format2 = format.upper()
131        if format2 == 'ASAP':
132            self._save(name)
133        else:
134            print "NYI"
135#             from asap._asap import sdwriter as _sw
136#             w = _sw(format2)
137#             w.write(self, name, stokes)
138
139        print_log()
140        return
141
142    def copy(self):
143        """
144        Return a copy of this scantable.
145        Parameters:
146            none
147        Example:
148            copiedscan = scan.copy()
149        """
150        sd = scantable(Scantable._copy(self))
151        return sd
152
153    def get_scan(self, scanid=None):
154        """
155        Return a specific scan (by scanno) or collection of scans (by
156        source name) in a new scantable.
157        Parameters:
158            scanid:    a (list of) scanno or a source name, unix-style
159                       patterns are accepted for source name matching, e.g.
160                       '*_R' gets all 'ref scans
161        Example:
162            # get all scans containing the source '323p459'
163            newscan = scan.get_scan('323p459')
164            # get all 'off' scans
165            refscans = scan.get_scan('*_R')
166            # get a susbset of scans by scanno (as listed in scan.summary())
167            newscan = scan.get_scan([0,2,7,10])
168        """
169        if scanid is None:
170            if rcParams['verbose']:
171                print "Please specify a scan no or name to retrieve from the scantable"
172                return
173            else:
174                raise RuntimeError("No scan given")
175
176        try:
177            bsel = self._getselection()
178            sel = asap._asap.Selector()
179            if type(scanid) is str:
180                sel._setname()
181                self._setselection(sel)
182                scopy = self._copy()
183                self._setselection(bsel)
184                return scantable(scopy)
185            elif type(scanid) is int:
186                sel._setscans([scanid])
187                self._setselection(sel)
188                scopy = self._copy()
189                self._setselection(bsel)
190                return scantable(scopy)
191            elif type(scanid) is list:
192                sel._setscans(scanid)
193                self._setselection(sel)
194                scopy = self._copy()
195                self._setselection(bsel)
196                return scantable(scopy)
197            else:
198                msg = "Illegal scanid type, use 'int' or 'list' if ints."
199                if rcParams['verbose']:
200                    print msg
201                else:
202                    raise TypeError(msg)
203        except RuntimeError:
204            if rcParams['verbose']: print "Couldn't find any match."
205            else: raise
206
207    def __str__(self):
208        return Scantable._summary(self,True)
209
210    def summary(self, filename=None, verbose=None):
211        """
212        Print a summary of the contents of this scantable.
213        Parameters:
214            filename:    the name of a file to write the putput to
215                         Default - no file output
216            verbose:     print extra info such as the frequency table
217                         The default (False) is taken from .asaprc
218        """
219        info = Scantable._summary(self, verbose)
220        if verbose is None: verbose = rcParams['scantable.verbosesummary']
221        if filename is not None:
222            if filename is "":
223                filename = 'scantable_summary.txt'
224            from os.path import expandvars, isdir
225            filename = expandvars(filename)
226            if not isdir(filename):
227                data = open(filename, 'w')
228                data.write(info)
229                data.close()
230            else:
231                msg = "Illegal file name '%s'." % (filename)
232                if rcParams['verbose']:
233                    print msg
234                else:
235                    raise IOError(msg)
236        if rcParams['verbose']:
237            try:
238                from IPython.genutils import page as pager
239            except ImportError:
240                from pydoc import pager
241            pager(info)
242        else:
243            return info
244
245    def set_cursor(self, beam=0, IF=0, pol=0):
246        """
247        Set the spectrum for individual operations.
248        Parameters:
249            beam, IF, pol:    a number
250        Example:
251            scan.set_cursor(0,0,1)
252            pol1sig = scan.stats(all=False) # returns std dev for beam=0
253                                            # if=0, pol=1
254        """
255        print "DEPRECATED"
256        varlist = vars()
257        sel = asap._asap.Selector()
258        sel._setbeams([beam])
259        sel._setpols([pol])
260        sel._setifs([IF])
261        self._add_history("set_cursor", varlist)
262        return
263
264    def get_cursor(self):
265        """
266        Return/print a the current 'cursor' into the Beam/IF/Pol cube.
267        Parameters:
268            none
269        Returns:
270            a list of values (currentBeam,currentIF,currentPol)
271        Example:
272            none
273        """
274        print "DEPRECATED"
275        sel = self._getselection()
276        i = sel.getbeams()[0]
277        j = sel.getifs()[0]
278        k = sel.getpols()[0]
279        from asap import asaplog
280        out = "--------------------------------------------------\n"
281        out += " Cursor position\n"
282        out += "--------------------------------------------------\n"
283        out += 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
284        asaplog.push(out)
285        print_log()
286        return i,j,k
287
288    def stats(self, stat='stddev', mask=None):
289        """
290        Determine the specified statistic of the current beam/if/pol
291        Takes a 'mask' as an optional parameter to specify which
292        channels should be excluded.
293        Parameters:
294            stat:    'min', 'max', 'sumsq', 'sum', 'mean'
295                     'var', 'stddev', 'avdev', 'rms', 'median'
296            mask:    an optional mask specifying where the statistic
297                     should be determined.
298        Example:
299            scan.set_unit('channel')
300            msk = scan.create_mask([100,200],[500,600])
301            scan.stats(stat='mean', mask=m)
302        """
303        from numarray import array,zeros,Float
304        if mask == None:
305            mask = []
306        axes = ['Beam','IF','Pol','Time']
307        if not self._check_ifs():
308             raise ValueError("Cannot apply mask as the IFs have different number of channels"
309                              "Please use setselection() to select individual IFs")
310
311        statvals = self._math._stats(self, mask, stat)
312        out = ''
313        axes = []
314        for i in range(self.nrow()):
315            axis = []
316            axis.append(self.getscan(i))
317            axis.append(self.getbeam(i))
318            axis.append(self.getif(i))
319            axis.append(self.getpol(i))
320            axis.append(self.getcycle(i))
321            axes.append(axis)
322            tm = self._gettime(i)
323            src = self._getsourcename(i)
324            out += 'Scan[%d] (%s) ' % (axis[0], src)
325            out += 'Time[%s]:\n' % (tm)
326            if self.nbeam(-1) > 1: out +=  ' Beam[%d] ' % (axis[1])
327            if self.nif(-1) > 1: out +=  ' IF[%d] ' % (axis[2])
328            if self.npol(-1) > 1: out +=  ' Pol[%d] ' % (axis[3])
329            out += '= %3.3f\n' % (statvals[i])
330            out +=  "--------------------------------------------------\n"
331
332        if rcParams['verbose']:
333            print "--------------------------------------------------"
334            print " ",stat
335            print "--------------------------------------------------"
336            print out
337        retval = { 'axesnames': ['scanno','beamno','ifno','polno','cycleno'],
338                   'axes' : axes,
339                   'data': statvals}
340        return retval
341
342    def stddev(self,mask=None):
343        """
344        Determine the standard deviation of the current beam/if/pol
345        Takes a 'mask' as an optional parameter to specify which
346        channels should be excluded.
347        Parameters:
348            mask:    an optional mask specifying where the standard
349                     deviation should be determined.
350
351        Example:
352            scan.set_unit('channel')
353            msk = scan.create_mask([100,200],[500,600])
354            scan.stddev(mask=m)
355        """
356        return self.stats(stat='stddev',mask=mask);
357
358    def get_tsys(self):
359        """
360        Return the System temperatures.
361        Parameters:
362
363        Returns:
364            a list of Tsys values for the current selection
365        """
366
367        return self._row_callback(self._gettsys, "Tsys")
368
369    def _row_callback(self, callback, label):
370        axes = []
371        axesnames = ['scanno','beamno','ifno','polno','cycleno']
372        out = ""
373        outvec =[]
374        for i in range(self.nrow()):
375            axis = []
376            axis.append(self.getscan(i))
377            axis.append(self.getbeam(i))
378            axis.append(self.getif(i))
379            axis.append(self.getpol(i))
380            axis.append(self.getcycle(i))
381            axes.append(axis)
382            tm = self._gettime(i)
383            src = self._getsourcename(i)
384            out += 'Scan[%d] (%s) ' % (axis[0], src)
385            out += 'Time[%s]:\n' % (tm)
386            if self.nbeam(-1) > 1: out +=  ' Beam[%d] ' % (axis[1])
387            if self.nif(-1) > 1: out +=  ' IF[%d] ' % (axis[2])
388            if self.npol(-1) > 1: out +=  ' Pol[%d] ' % (axis[3])
389            outvec.append(callback(i))
390            out += '= %3.3f\n' % (outvec[i])
391            out +=  "--------------------------------------------------\n"
392        if rcParams['verbose']:
393            print "--------------------------------------------------"
394            print " %s" % (label)
395            print "--------------------------------------------------"
396            print out
397        retval = {'axesnames': axesnames, 'axes': axes, 'data': outvec}
398        return retval
399
400
401    def get_time(self, row=-1):
402        """
403        Get a list of time stamps for the observations.
404        Return a string for each integration in the scantable.
405        Parameters:
406            row:    row no of integration. Default -1 return all rows
407        Example:
408            none
409        """
410        out = []
411        if row == -1:
412            for i in range(self.nrow()):
413                out.append(self._gettime(i))
414            return out
415        else:
416            if row < self.nrow():
417                return self._gettime(row)
418
419    def get_sourcename(self, row=-1):
420        """
421        Get a list source names for the observations.
422        Return a string for each integration in the scantable.
423        Parameters:
424            row:    row no of integration. Default -1 return all rows
425        Example:
426            none
427        """
428        out = []
429        if row == -1:
430            return [self._getsourcename(i) for i in range(self.nrow())]
431        else:
432            if  0 <= row < self.nrow():
433                return self._getsourcename(row)
434
435    def get_elevation(self, row=-1):
436        """
437        Get a list of elevations for the observations.
438        Return a float for each integration in the scantable.
439        Parameters:
440            row:    row no of integration. Default -1 return all rows
441        Example:
442            none
443        """
444        out = []
445        if row == -1:
446            return [self._getelevation(i) for i in range(self.nrow())]
447        else:
448            if  0 <= row < self.nrow():
449                return self._getelevation(row)
450
451    def get_azimuth(self, row=-1):
452        """
453        Get a list of azimuths for the observations.
454        Return a float for each integration in the scantable.
455        Parameters:
456            row:    row no of integration. Default -1 return all rows
457        Example:
458            none
459        """
460        out = []
461        if row == -1:
462            return [self._getazimuth(i) for i in range(self.nrow())]
463        else:
464            if  0 <= row < self.nrow():
465                return self._getazimuth(row)
466
467    def get_parangle(self, row=-1):
468        """
469        Get a list of parallactic angles for the observations.
470        Return a float for each integration in the scantable.
471        Parameters:
472            row:    row no of integration. Default -1 return all rows
473        Example:
474            none
475        """
476        out = []
477        if row == -1:
478            return [self._getparangle(i) for i in range(self.nrow())]
479        else:
480            if  0 <= row < self.nrow():
481                return self._getparangle(row)
482
483    def set_unit(self, unit='channel'):
484        """
485        Set the unit for all following operations on this scantable
486        Parameters:
487            unit:    optional unit, default is 'channel'
488                     one of '*Hz','km/s','channel', ''
489        """
490        varlist = vars()
491        if unit in ['','pixel', 'channel']:
492            unit = ''
493        inf = list(self._getcoordinfo())
494        inf[0] = unit
495        self._setcoordinfo(inf)
496        self._add_history("set_unit",varlist)
497
498    def set_instrument(self, instr):
499        """
500        Set the instrument for subsequent processing
501        Parameters:
502            instr:    Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
503                      'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
504        """
505        self._setInstrument(instr)
506        self._add_history("set_instument",vars())
507        print_log()
508
509    def set_doppler(self, doppler='RADIO'):
510        """
511        Set the doppler for all following operations on this scantable.
512        Parameters:
513            doppler:    One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
514        """
515        varlist = vars()
516        inf = list(self._getcoordinfo())
517        inf[2] = doppler
518        self._setcoordinfo(inf)
519        self._add_history("set_doppler",vars())
520        print_log()
521
522    def set_freqframe(self, frame=None):
523        """
524        Set the frame type of the Spectral Axis.
525        Parameters:
526            frame:   an optional frame type, default 'LSRK'. Valid frames are:
527                     'REST','TOPO','LSRD','LSRK','BARY',
528                     'GEO','GALACTO','LGROUP','CMB'
529        Examples:
530            scan.set_freqframe('BARY')
531        """
532        if frame is None: frame = rcParams['scantable.freqframe']
533        varlist = vars()
534        valid = ['REST','TOPO','LSRD','LSRK','BARY', \
535                   'GEO','GALACTO','LGROUP','CMB']
536
537        if 1:#frame in valid:
538            inf = list(self._getcoordinfo())
539            inf[1] = frame
540            self._setcoordinfo(inf)
541            self._add_history("set_freqframe",varlist)
542        else:
543            msg  = "Please specify a valid freq type. Valid types are:\n",valid
544            if rcParams['verbose']:
545                print msg
546            else:
547                raise TypeError(msg)
548        print_log()
549
550    def get_unit(self):
551        """
552        Get the default unit set in this scantable
553        Parameters:
554        Returns:
555            A unit string
556        """
557        inf = self._getcoordinfo()
558        unit = inf[0]
559        if unit == '': unit = 'channel'
560        return unit
561
562    def get_abcissa(self, rowno=0):
563        """
564        Get the abcissa in the current coordinate setup for the currently
565        selected Beam/IF/Pol
566        Parameters:
567            rowno:    an optional row number in the scantable. Default is the
568                      first row, i.e. rowno=0
569        Returns:
570            The abcissa values and it's format string (as a dictionary)
571        """
572        abc = self._getabcissa(rowno)
573        lbl = self._getabcissalabel(rowno)
574        print_log()
575        return abc, lbl
576
577    def create_mask(self, *args, **kwargs):
578        """
579        Compute and return a mask based on [min,max] windows.
580        The specified windows are to be INCLUDED, when the mask is
581        applied.
582        Parameters:
583            [min,max],[min2,max2],...
584                Pairs of start/end points specifying the regions
585                to be masked
586            invert:     optional argument. If specified as True,
587                        return an inverted mask, i.e. the regions
588                        specified are EXCLUDED
589            row:        create the mask using the specified row for
590                        unit conversions, default is row=0
591                        only necessary if frequency varies over rows.
592        Example:
593            scan.set_unit('channel')
594
595            a)
596            msk = scan.set_mask([400,500],[800,900])
597            # masks everything outside 400 and 500
598            # and 800 and 900 in the unit 'channel'
599
600            b)
601            msk = scan.set_mask([400,500],[800,900], invert=True)
602            # masks the regions between 400 and 500
603            # and 800 and 900 in the unit 'channel'
604
605        """
606        row = 0
607        if kwargs.has_key("row"):
608            row = kwargs.get("row")
609        data = self._getabcissa(row)
610        u = self._getcoordinfo()[0]
611        if rcParams['verbose']:
612            if u == "": u = "channel"
613            from asap import asaplog
614            msg = "The current mask window unit is %s" % u
615            if not self._check_ifs():
616                msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
617            asaplog.push(msg)
618        n = self.nchan()
619        msk = zeros(n)
620        # test if args is a 'list' or a 'normal *args - UGLY!!!
621
622        ws = (isinstance(args[-1][-1],int) or isinstance(args[-1][-1],float)) and args or args[0]
623        for window in ws:
624            if (len(window) != 2 or window[0] > window[1] ):
625                raise TypeError("A window needs to be defined as [min,max]")
626            for i in range(n):
627                if data[i] >= window[0] and data[i] < window[1]:
628                    msk[i] = 1
629        if kwargs.has_key('invert'):
630            if kwargs.get('invert'):
631                from numarray import logical_not
632                msk = logical_not(msk)
633        print_log()
634        return msk
635
636    def get_restfreqs(self):
637        """
638        Get the restfrequency(s) stored in this scantable.
639        The return value(s) are always of unit 'Hz'
640        Parameters:
641            none
642        Returns:
643            a list of doubles
644        """
645        return list(self._getrestfreqs())
646
647
648    def set_restfreqs(self, freqs=None, unit='Hz'):
649        """
650        Set or replace the restfrequency specified and
651        If the 'freqs' argument holds a scalar,
652        then that rest frequency will be applied to all the selected
653        data.  If the 'freqs' argument holds
654        a vector, then it MUST be of equal or smaller length than
655        the number of IFs (and the available restfrequencies will be
656        replaced by this vector).  In this case, *all* data have
657        the restfrequency set per IF according
658        to the corresponding value you give in the 'freqs' vector.
659        E.g. 'freqs=[1e9,2e9]'  would mean IF 0 gets restfreq 1e9 and
660        IF 1 gets restfreq 2e9.
661        You can also specify the frequencies via known line names
662        from the built-in Lovas table.
663        Parameters:
664            freqs:   list of rest frequency values or string idenitfiers
665            unit:    unit for rest frequency (default 'Hz')
666
667        Example:
668            # set the given restfrequency for the whole table
669            scan.set_restfreqs(freqs=1.4e9)
670            # If thee number of IFs in the data is >= 2 the IF0 gets the first
671            # value IF1 the second...
672            scan.set_restfreqs(freqs=[1.4e9,1.67e9])
673            #set the given restfrequency for the whole table (by name)
674            scan.set_restfreqs(freqs="OH1667")
675
676        Note:
677            To do more sophisticate Restfrequency setting, e.g. on a
678            source and IF basis, use scantable.set_selection() before using
679            this function.
680            # provide your scantable is call scan
681            selection = selector()
682            selection.set_name("ORION*")
683            selection.set_ifs([1])
684            scan.set_selection(selection)
685            scan.set_restfreqs(freqs=86.6e9)
686
687        """
688        varlist = vars()
689
690        t = type(freqs)
691        if isinstance(freqs, int) or isinstance(freqs,float):
692           self._setrestfreqs(freqs, unit)
693        elif isinstance(freqs, list) or isinstance(freqs,tuple):
694            if isinstance(freqs[-1], int) or isinstance(freqs[-1],float):
695                from asap._asap import selector
696                sel = selector()
697                savesel = self._getselection()
698                for i in xrange(len(freqs)):
699                    sel._setifs([i])
700                    self._setselection(sel)
701                    self._setrestfreqs(freqs[i], unit)
702                self._setselection(savesel)
703            elif isinstance(freqs[-1], str):
704                # not yet implemented
705                pass
706        else:
707            return
708        self._add_history("set_restfreqs", varlist)
709
710
711
712    def history(self):
713        hist = list(self._gethistory())
714        out = "-"*80
715        for h in hist:
716            if h.startswith("---"):
717                out += "\n"+h
718            else:
719                items = h.split("##")
720                date = items[0]
721                func = items[1]
722                items = items[2:]
723                out += "\n"+date+"\n"
724                out += "Function: %s\n  Parameters:" % (func)
725                for i in items:
726                    s = i.split("=")
727                    out += "\n   %s = %s" % (s[0],s[1])
728                out += "\n"+"-"*80
729        try:
730            from IPython.genutils import page as pager
731        except ImportError:
732            from pydoc import pager
733        pager(out)
734        return
735
736    #
737    # Maths business
738    #
739
740    def average_time(self, mask=None, scanav=False, weight='tint', align=False):
741        """
742        Return the (time) average of a scan, or apply it 'insitu'.
743        Note:
744            in channels only
745            The cursor of the output scan is set to 0.
746        Parameters:
747            one scan or comma separated  scans
748            mask:     an optional mask (only used for 'var' and 'tsys'
749                      weighting)
750            scanav:   True averages each scan separately
751                      False (default) averages all scans together,
752            weight:   Weighting scheme. 'none', 'var' (1/var(spec)
753                      weighted), 'tsys' (1/Tsys**2 weighted), 'tint'
754                      (integration time weighted) or 'tintsys' (Tint/Tsys**2).
755                      The default is 'tint'
756            align:    align the spectra in velocity before averaging. It takes
757                      the time of the first spectrum as reference time.
758        Example:
759            # time average the scantable without using a mask
760            newscan = scan.average_time()
761        """
762        varlist = vars()
763        if weight is None: weight = 'tint'
764        if mask is None: mask = ()
765        if scanav:
766          scanav = "SCAN"
767        else:
768          scanav = "NONE"
769        s = scantable(self._math._average((self,), mask, weight.upper(),
770                      scanav, align))
771        s._add_history("average_time",varlist)
772        print_log()
773        return s
774
775    def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
776        """
777        Return a scan where all spectra are converted to either
778        Jansky or Kelvin depending upon the flux units of the scan table.
779        By default the function tries to look the values up internally.
780        If it can't find them (or if you want to over-ride), you must
781        specify EITHER jyperk OR eta (and D which it will try to look up
782        also if you don't set it). jyperk takes precedence if you set both.
783        Parameters:
784            jyperk:      the Jy / K conversion factor
785            eta:         the aperture efficiency
786            d:           the geomtric diameter (metres)
787            insitu:      if False a new scantable is returned.
788                         Otherwise, the scaling is done in-situ
789                         The default is taken from .asaprc (False)
790            allaxes:         if True apply to all spectra. Otherwise
791                         apply only to the selected (beam/pol/if)spectra only
792                         The default is taken from .asaprc (True if none)
793        """
794        if insitu is None: insitu = rcParams['insitu']
795        self._math._setinsitu(insitu)
796        varlist = vars()
797        if jyperk is None: jyperk = -1.0
798        if d is None: d = -1.0
799        if eta is None: eta = -1.0
800        s = scantable(self._math._convertflux(self, d, eta, jyperk))
801        s._add_history("convert_flux", varlist)
802        print_log()
803        if insitu: self._assign(s)
804        else: return s
805
806    def gain_el(self, poly=None, filename="", method="linear", insitu=None):
807        """
808        Return a scan after applying a gain-elevation correction.
809        The correction can be made via either a polynomial or a
810        table-based interpolation (and extrapolation if necessary).
811        You specify polynomial coefficients, an ascii table or neither.
812        If you specify neither, then a polynomial correction will be made
813        with built in coefficients known for certain telescopes (an error
814        will occur if the instrument is not known).
815        The data and Tsys are *divided* by the scaling factors.
816        Parameters:
817            poly:        Polynomial coefficients (default None) to compute a
818                         gain-elevation correction as a function of
819                         elevation (in degrees).
820            filename:    The name of an ascii file holding correction factors.
821                         The first row of the ascii file must give the column
822                         names and these MUST include columns
823                         "ELEVATION" (degrees) and "FACTOR" (multiply data
824                         by this) somewhere.
825                         The second row must give the data type of the
826                         column. Use 'R' for Real and 'I' for Integer.
827                         An example file would be
828                         (actual factors are arbitrary) :
829
830                         TIME ELEVATION FACTOR
831                         R R R
832                         0.1 0 0.8
833                         0.2 20 0.85
834                         0.3 40 0.9
835                         0.4 60 0.85
836                         0.5 80 0.8
837                         0.6 90 0.75
838            method:      Interpolation method when correcting from a table.
839                         Values are  "nearest", "linear" (default), "cubic"
840                         and "spline"
841            insitu:      if False a new scantable is returned.
842                         Otherwise, the scaling is done in-situ
843                         The default is taken from .asaprc (False)
844        """
845
846        if insitu is None: insitu = rcParams['insitu']
847        self._math._setinsitu(insitu)
848        varlist = vars()
849        if poly is None:
850           poly = ()
851        from os.path import expandvars
852        filename = expandvars(filename)
853        s = scantable(self._math._gainel(self, poly, filename, method))
854        s._add_history("gain_el", varlist)
855        print_log()
856        if insitu: self._assign(s)
857        else: return s
858
859    def freq_align(self, reftime=None, method='cubic', insitu=None):
860        """
861        Return a scan where all rows have been aligned in frequency/velocity.
862        The alignment frequency frame (e.g. LSRK) is that set by function
863        set_freqframe.
864        Parameters:
865            reftime:     reference time to align at. By default, the time of
866                         the first row of data is used.
867            method:      Interpolation method for regridding the spectra.
868                         Choose from "nearest", "linear", "cubic" (default)
869                         and "spline"
870            insitu:      if False a new scantable is returned.
871                         Otherwise, the scaling is done in-situ
872                         The default is taken from .asaprc (False)
873        """
874        if insitu is None: insitu = rcParams["insitu"]
875        self._math._setinsitu(insitu)
876        varlist = vars()
877        if reftime is None: reftime = ""
878        s = scantable(self._math._freq_align(self, reftime, method))
879        s._add_history("freq_align", varlist)
880        print_log()
881        if insitu: self._assign(s)
882        else: return s
883
884    def opacity(self, tau, insitu=None):
885        """
886        Apply an opacity correction. The data
887        and Tsys are multiplied by the correction factor.
888        Parameters:
889            tau:         Opacity from which the correction factor is
890                         exp(tau*ZD)
891                         where ZD is the zenith-distance
892            insitu:      if False a new scantable is returned.
893                         Otherwise, the scaling is done in-situ
894                         The default is taken from .asaprc (False)
895        """
896        if insitu is None: insitu = rcParams['insitu']
897        self._math._setinsitu(insitu)
898        varlist = vars()
899        s = scantable(self._math._opacity(self, tau))
900        s._add_history("opacity", varlist)
901        print_log()
902        if insitu: self._assign(s)
903        else: return s
904
905    def bin(self, width=5, insitu=None):
906        """
907        Return a scan where all spectra have been binned up.
908            width:       The bin width (default=5) in pixels
909            insitu:      if False a new scantable is returned.
910                         Otherwise, the scaling is done in-situ
911                         The default is taken from .asaprc (False)
912        """
913        if insitu is None: insitu = rcParams['insitu']
914        self._math._setinsitu(insitu)
915        varlist = vars()
916        s = scantable(self._math._bin(self, width))
917        s._add_history("bin",varlist)
918        print_log()
919        if insitu: self._assign(s)
920        else: return s
921
922
923    def resample(self, width=5, method='cubic', insitu=None):
924        """
925        Return a scan where all spectra have been binned up
926            width:       The bin width (default=5) in pixels
927            method:      Interpolation method when correcting from a table.
928                         Values are  "nearest", "linear", "cubic" (default)
929                         and "spline"
930            insitu:      if False a new scantable is returned.
931                         Otherwise, the scaling is done in-situ
932                         The default is taken from .asaprc (False)
933        """
934        if insitu is None: insitu = rcParams['insitu']
935        self._math._setinsitu(insitu)
936        varlist = vars()
937        s = scantable(self._math._resample(self, method, width))
938        s._add_history("resample",varlist)
939        print_log()
940        if insitu: self._assign(s)
941        else: return s
942
943
944#     def average_pol(self, mask=None, weight='none', insitu=None):
945#         """
946#         Average the Polarisations together.
947#         The polarisation cursor of the output scan is set to 0
948#         Parameters:
949#             mask:        An optional mask defining the region, where the
950#                          averaging will be applied. The output will have all
951#                          specified points masked.
952#             weight:      Weighting scheme. 'none' (default), 'var' (1/var(spec)
953#                          weighted), or 'tsys' (1/Tsys**2 weighted)
954#             insitu:      if False a new scantable is returned.
955#                          Otherwise, the scaling is done in-situ
956#                          The default is taken from .asaprc (False)
957#         """
958#         if insitu is None: insitu = rcParams['insitu']
959#         self._math._setinsitu(insitu)
960#         varlist = vars()
961#
962#         if mask is None:
963#         mask = ()
964#         s = self._math._averagepol(self, mask, weight)
965#         s._add_history("average_pol",varlist)
966#         print_log()
967#         if insitu: self._assign(s)
968#         else: return scantable(s)
969
970    def smooth(self, kernel="hanning", width=5.0, insitu=None):
971        """
972        Smooth the spectrum by the specified kernel (conserving flux).
973        Parameters:
974            scan:       The input scan
975            kernel:     The type of smoothing kernel. Select from
976                        'hanning' (default), 'gaussian' and 'boxcar'.
977                        The first three characters are sufficient.
978            width:      The width of the kernel in pixels. For hanning this is
979                        ignored otherwise it defauls to 5 pixels.
980                        For 'gaussian' it is the Full Width Half
981                        Maximum. For 'boxcar' it is the full width.
982            insitu:     if False a new scantable is returned.
983                        Otherwise, the scaling is done in-situ
984                        The default is taken from .asaprc (False)
985        Example:
986             none
987        """
988        if insitu is None: insitu = rcParams['insitu']
989        self._math._setinsitu(insitu)
990        varlist = vars()
991        s = scantable(self._math._smooth(self,kernel,width))
992        s._add_history("smooth", varlist)
993        print_log()
994        if insitu: self._assign(s)
995        else: return s
996
997
998    def poly_baseline(self, mask=None, order=0, insitu=None):
999        """
1000        Return a scan which has been baselined (all rows) by a polynomial.
1001        Parameters:
1002            scan:       a scantable
1003            mask:       an optional mask
1004            order:      the order of the polynomial (default is 0)
1005            insitu:     if False a new scantable is returned.
1006                        Otherwise, the scaling is done in-situ
1007                        The default is taken from .asaprc (False)
1008            allaxes:    If True (default) apply to all spectra. Otherwise
1009                        apply only to the selected (beam/pol/if)spectra only
1010                        The default is taken from .asaprc (True if none)
1011        Example:
1012            # return a scan baselined by a third order polynomial,
1013            # not using a mask
1014            bscan = scan.poly_baseline(order=3)
1015        """
1016        if insitu is None: insitu = rcParams['insitu']
1017        varlist = vars()
1018        if mask is None:
1019            from numarray import ones
1020            mask = list(ones(self.nchan(-1)))
1021        from asap.asapfitter import fitter
1022        f = fitter()
1023        f.set_scan(self, mask)
1024        f.set_function(poly=order)
1025        s = f.auto_fit(insitu)
1026        s._add_history("poly_baseline", varlist)
1027        print_log()
1028        if insitu: self._assign(s)
1029        else: return s
1030
1031    def auto_poly_baseline(self, mask=[], edge=(0,0), order=0,
1032                           threshold=3, insitu=None):
1033        """
1034        Return a scan which has been baselined (all rows) by a polynomial.
1035        Spectral lines are detected first using linefinder and masked out
1036        to avoid them affecting the baseline solution.
1037
1038        Parameters:
1039            mask:       an optional mask retreived from scantable
1040            edge:       an optional number of channel to drop at
1041                        the edge of spectrum. If only one value is
1042                        specified, the same number will be dropped from
1043                        both sides of the spectrum. Default is to keep
1044                        all channels. Nested tuples represent individual
1045                        edge selection for different IFs (a number of spectral
1046                        channels can be different)
1047            order:      the order of the polynomial (default is 0)
1048            threshold:  the threshold used by line finder. It is better to
1049                        keep it large as only strong lines affect the
1050                        baseline solution.
1051            insitu:     if False a new scantable is returned.
1052                        Otherwise, the scaling is done in-situ
1053                        The default is taken from .asaprc (False)
1054
1055        Example:
1056            scan2=scan.auto_poly_baseline(order=7)
1057        """
1058        if insitu is None: insitu = rcParams['insitu']
1059        varlist = vars()
1060        from asap.asapfitter import fitter
1061        from asap.asaplinefind import linefinder
1062        from asap import _is_sequence_or_number as _is_valid
1063
1064        # check whether edge is set up for each IF individually
1065        individualEdge = False;
1066        if len(edge)>1:
1067           if isinstance(edge[0],list) or isinstance(edge[0],tuple):
1068               individualEdge = True;
1069
1070        if not _is_valid(edge, int) and not individualEdge:
1071            raise ValueError, "Parameter 'edge' has to be an integer or a \
1072            pair of integers specified as a tuple. Nested tuples are allowed \
1073            to make individual selection for different IFs."
1074
1075        curedge = (0,0)
1076        if individualEdge:
1077           for edge_par in edge:
1078               if not _is_valid(edge,int):
1079                  raise ValueError, "Each element of the 'edge' tuple has \
1080                  to be a pair of integers or an integer."
1081        else:
1082           curedge = edge;
1083
1084        # setup fitter
1085        f = fitter()
1086        f.set_function(poly=order)
1087
1088        # setup line finder
1089        fl=linefinder()
1090        fl.set_options(threshold=threshold)
1091
1092        if not insitu:
1093            workscan=self.copy()
1094        else:
1095            workscan=self
1096
1097        fl.set_scan(workscan)
1098
1099        rows=range(workscan.nrow())
1100        from asap import asaplog
1101        asaplog.push("Processing:")
1102        for r in rows:
1103            msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" %        (workscan.getscan(r),workscan.getbeam(r),workscan.getif(r),workscan.getpol(r), workscan.getcycle(r))
1104            asaplog.push(msg, False)
1105
1106            # figure out edge parameter
1107            if individualEdge:
1108               if len(edge)>=workscan.getif(r):
1109                  raise RuntimeError, "Number of edge elements appear to be less than the number of IFs"
1110                  curedge = edge[workscan.getif(r)]
1111
1112            # setup line finder
1113            fl.find_lines(r,mask,curedge)
1114            f.set_scan(workscan, fl.get_mask())
1115            f.x = workscan._getabcissa(r)
1116            f.y = workscan._getspectrum(r)
1117            f.data = None
1118            f.fit()
1119            x = f.get_parameters()
1120            workscan._setspectrum(f.fitter.getresidual(), r)
1121        workscan._add_history("poly_baseline", varlist)
1122        if insitu:
1123            self._assign(workscan)
1124        else:
1125            return workscan
1126
1127    def rotate_linpolphase(self, angle):
1128        """
1129        Rotate the phase of the complex polarization O=Q+iU correlation.
1130        This is always done in situ in the raw data.  So if you call this
1131        function more than once then each call rotates the phase further.
1132        Parameters:
1133            angle:   The angle (degrees) to rotate (add) by.
1134        Examples:
1135            scan.rotate_linpolphase(2.3)
1136        """
1137        varlist = vars()
1138        self._math._rotate_linpolphase(self, angle)
1139        self._add_history("rotate_linpolphase", varlist)
1140        print_log()
1141        return
1142
1143
1144    def rotate_xyphase(self, angle):
1145        """
1146        Rotate the phase of the XY correlation.  This is always done in situ
1147        in the data.  So if you call this function more than once
1148        then each call rotates the phase further.
1149        Parameters:
1150            angle:   The angle (degrees) to rotate (add) by.
1151        Examples:
1152            scan.rotate_xyphase(2.3)
1153        """
1154        varlist = vars()
1155        self._math._rotate_xyphase(self, angle)
1156        self._add_history("rotate_xyphase", varlist)
1157        print_log()
1158        return
1159
1160    def swap_linears(self):
1161        """
1162        Swap the linear polarisations XX and YY
1163        """
1164        varlist = vars()
1165        self._math._swap_linears(self)
1166        self._add_history("swap_linears", varlist)
1167        print_log()
1168        return
1169
1170    def invert_phase(self):
1171        """
1172        Invert the phase of the complex polarisation
1173        """
1174        varlist = vars()
1175        self._math._invert_phase(self)
1176        self._add_history("invert_phase", varlist)
1177        print_log()
1178        return
1179
1180    def add(self, offset, insitu=None):
1181        """
1182        Return a scan where all spectra have the offset added
1183        Parameters:
1184            offset:      the offset
1185            insitu:      if False a new scantable is returned.
1186                         Otherwise, the scaling is done in-situ
1187                         The default is taken from .asaprc (False)
1188        """
1189        if insitu is None: insitu = rcParams['insitu']
1190        self._math._setinsitu(insitu)
1191        varlist = vars()
1192        s = scantable(self._math._unaryop(self, offset, "ADD", False))
1193        s._add_history("add",varlist)
1194        print_log()
1195        if insitu:
1196            self._assign(s)
1197        else:
1198            return s
1199
1200    def scale(self, factor, tsys=True, insitu=None,):
1201        """
1202        Return a scan where all spectra are scaled by the give 'factor'
1203        Parameters:
1204            factor:      the scaling factor
1205            insitu:      if False a new scantable is returned.
1206                         Otherwise, the scaling is done in-situ
1207                         The default is taken from .asaprc (False)
1208            tsys:        if True (default) then apply the operation to Tsys
1209                         as well as the data
1210        """
1211        if insitu is None: insitu = rcParams['insitu']
1212        self._math._setinsitu(insitu)
1213        varlist = vars()
1214        s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
1215        s._add_history("scale",varlist)
1216        print_log()
1217        if insitu:
1218            self._assign(s)
1219        else:
1220            return s
1221
1222    def auto_quotient(self, mode='time', preserve=True):
1223        """
1224        This function allows to build quotients automatically.
1225        It assumes the observation to have the same numer of
1226        "ons" and "offs"
1227        It will support "closest off in time" in the future
1228        Parameters:
1229            mode:           the on/off detection mode; 'suffix' (default)
1230                            'suffix' identifies 'off' scans by the
1231                            trailing '_R' (Mopra/Parkes) or
1232                            '_e'/'_w' (Tid)
1233            preserve:       you can preserve (default) the continuum or
1234                            remove it.  The equations used are
1235                            preserve: Output = Toff * (on/off) - Toff
1236                            remove:   Output = Tref * (on/off) - Ton
1237        """
1238        modes = ["time"]
1239        if not mode in modes:
1240            msg = "please provide valid mode. Valid modes are %s" % (modes)
1241            raise ValueError(msg)
1242        varlist = vars()
1243        s = scantable(self._math._quotient(self, mode, preserve))
1244        s._add_history("auto_quotient",varlist)
1245        print_log()
1246        return s
1247
1248
1249
1250
1251    def freq_switch(self, insitu=None):
1252        """
1253        Apply frequency switching to the data.
1254        Parameters:
1255            insitu:      if False a new scantable is returned.
1256                         Otherwise, the swictching is done in-situ
1257                         The default is taken from .asaprc (False)
1258        Example:
1259            none
1260        """
1261        if insitu is None: insitu = rcParams['insitu']
1262        self._math._setinsitu(insitu)
1263        varlist = vars()
1264        s = scantable(self._math._freqswitch(self))
1265        s._add_history("freq_switch",varlist)
1266        print_log()
1267        if insitu: self._assign(s)
1268        else: return s
1269
1270    def recalc_azel(self):
1271        """
1272        Recalculate the azimuth and elevation for each position.
1273        Parameters:
1274            none
1275        Example:
1276        """
1277        varlist = vars()
1278        self._recalcazel()
1279        self._add_history("recalc_azel", varlist)
1280        print_log()
1281        return
1282
1283    def __add__(self, other):
1284        varlist = vars()
1285        s = None
1286        if isinstance(other, scantable):
1287            print "scantable + scantable NYI"
1288            return
1289        elif isinstance(other, float):
1290            s = scantable(self._math._unaryop(self, other, "ADD", False))
1291        else:
1292            raise TypeError("Other input is not a scantable or float value")
1293        s._add_history("operator +", varlist)
1294        print_log()
1295        return s
1296
1297    def __sub__(self, other):
1298        """
1299        implicit on all axes and on Tsys
1300        """
1301        varlist = vars()
1302        s = None
1303        if isinstance(other, scantable):
1304            print "scantable - scantable NYI"
1305            return
1306        elif isinstance(other, float):
1307            s = scantable(self._math._unaryop(self, other, "SUB", False))
1308        else:
1309            raise TypeError("Other input is not a scantable or float value")
1310        s._add_history("operator -", varlist)
1311        print_log()
1312        return s
1313
1314    def __mul__(self, other):
1315        """
1316        implicit on all axes and on Tsys
1317        """
1318        varlist = vars()
1319        s = None
1320        if isinstance(other, scantable):
1321            print "scantable * scantable NYI"
1322            return
1323        elif isinstance(other, float):
1324            s = scantable(self._math._unaryop(self, other, "MUL", False))
1325        else:
1326            raise TypeError("Other input is not a scantable or float value")
1327        s._add_history("operator *", varlist)
1328        print_log()
1329        return s
1330
1331
1332    def __div__(self, other):
1333        """
1334        implicit on all axes and on Tsys
1335        """
1336        varlist = vars()
1337        s = None
1338        if isinstance(other, scantable):
1339            print "scantable / scantable NYI"
1340            return
1341        elif isinstance(other, float):
1342            if other == 0.0:
1343                raise ZeroDivisionError("Dividing by zero is not recommended")
1344            s = scantable(self._math._unaryop(self, other, "DIV", False))
1345        else:
1346            raise TypeError("Other input is not a scantable or float value")
1347        s._add_history("operator /", varlist)
1348        print_log()
1349        return s
1350
1351    def get_fit(self, row=0):
1352        """
1353        Print or return the stored fits for a row in the scantable
1354        Parameters:
1355            row:    the row which the fit has been applied to.
1356        """
1357        if row > self.nrow():
1358            return
1359        from asap import asapfit
1360        fit = asapfit(self._getfit(row))
1361        if rcParams['verbose']:
1362            print fit
1363            return
1364        else:
1365            return fit.as_dict()
1366
1367    def _add_history(self, funcname, parameters):
1368        # create date
1369        sep = "##"
1370        from datetime import datetime
1371        dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1372        hist = dstr+sep
1373        hist += funcname+sep#cdate+sep
1374        if parameters.has_key('self'): del parameters['self']
1375        for k,v in parameters.iteritems():
1376            if type(v) is dict:
1377                for k2,v2 in v.iteritems():
1378                    hist += k2
1379                    hist += "="
1380                    if isinstance(v2,scantable):
1381                        hist += 'scantable'
1382                    elif k2 == 'mask':
1383                        if isinstance(v2,list) or isinstance(v2,tuple):
1384                            hist += str(self._zip_mask(v2))
1385                        else:
1386                            hist += str(v2)
1387                    else:
1388                        hist += str(v2)
1389            else:
1390                hist += k
1391                hist += "="
1392                if isinstance(v,scantable):
1393                    hist += 'scantable'
1394                elif k == 'mask':
1395                    if isinstance(v,list) or isinstance(v,tuple):
1396                        hist += str(self._zip_mask(v))
1397                    else:
1398                        hist += str(v)
1399                else:
1400                    hist += str(v)
1401            hist += sep
1402        hist = hist[:-2] # remove trailing '##'
1403        self._addhistory(hist)
1404
1405
1406    def _zip_mask(self, mask):
1407        mask = list(mask)
1408        i = 0
1409        segments = []
1410        while mask[i:].count(1):
1411            i += mask[i:].index(1)
1412            if mask[i:].count(0):
1413                j = i + mask[i:].index(0)
1414            else:
1415                j = len(mask)
1416            segments.append([i,j])
1417            i = j
1418        return segments
1419
1420    def _get_ordinate_label(self):
1421        fu = "("+self.get_fluxunit()+")"
1422        import re
1423        lbl = "Intensity"
1424        if re.match(".K.",fu):
1425            lbl = "Brightness Temperature "+ fu
1426        elif re.match(".Jy.",fu):
1427            lbl = "Flux density "+ fu
1428        return lbl
1429
1430    def _check_ifs(self):
1431        nchans = [self.nchan(i) for i in range(self.nif(-1))]
1432        nchans = filter(lambda t: t > 0, nchans)
1433        return (sum(nchans)/len(nchans) == nchans[0])
Note: See TracBrowser for help on using the repository browser.