source: trunk/python/scantable.py @ 981

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

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