source: trunk/python/scantable.py @ 196

Last change on this file since 196 was 196, checked in by kil064, 19 years ago

imporve docs for function 'save'

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.9 KB
Line 
1from asap._asap import sdtable
2from numarray import ones,zeros
3import sys
4
5class scantable(sdtable):
6    """
7        The ASAP container for scans
8    """
9   
10    def __init__(self, filename):
11        """
12        Create a scantable from a saved one or make a reference
13        Parameters:
14            filename:    the name of an asap table on disk
15                         or
16                         the name of a rpfits/sdfits/ms file
17                         (integrations within scans are auto averaged
18                         and the whole file is read)
19                         or
20                         [advanced] a reference to an existing
21                         scantable
22        """
23        self._vb = True
24        self._p = None
25        from os import stat as st
26        import stat
27        if isinstance(filename,sdtable):
28            sdtable.__init__(self, filename)           
29        else:
30            mode = st(filename)[stat.ST_MODE]
31            if stat.S_ISDIR(mode):
32                # crude check if asap table
33                if stat.S_ISREG(st(filename+'/table.info')[stat.ST_MODE]):
34                    sdtable.__init__(self, filename)
35                else:
36                    print 'The given file is not a valid asap table'
37            else:           
38                from asap._asap import sdreader
39                r = sdreader(filename)
40                print 'Importing data...'
41                r.read([-1])
42                tbl = r.getdata()
43               
44                from asap._asap import average
45                tmp = tuple([tbl])
46                print 'Auto averging integrations...'
47                tbl2 = average(tmp,(),True,'none')
48                sdtable.__init__(self,tbl2)
49                del r,tbl
50
51    def save(self, name="", format='ASAP'):
52        """
53        Store the scantable on disk. This can be a asap file or SDFITS/MS2.
54        Parameters:
55            name:        the name of the outputfile. For format="FITS" this
56                         is the root file name (and can be empty).
57            format:      an optional file format. Default is ASAP.
58                         Allowed are - 'ASAP' (save as ASAP Table),
59                                       'SDFITS' (save as SDFITS file)
60                                       'FITS' (saves each row as a FITS Image)
61                                       'MS2' (saves as an aips++ MeasurementSet V2)
62        Example:
63            scan.save('myscan.asap')
64            scan.save('myscan.sdfits','SDFITS')
65        """
66        if format == 'ASAP':
67            self._save(name)
68        else:
69            from asap._asap import sdwriter as _sw
70            w = _sw(format)
71            w.write(self, name)
72        return
73
74    def _verbose(self, *args):
75        """
76        Set the verbose level True or False, to indicate if output
77        should be printed as well as returned.
78        """
79        if len(args) == 0:
80            return self._vb
81        elif type(args[0]) is bool:
82            self._vb = args[0]
83            return
84
85    def copy(self):
86        """
87        Return a copy of this scantable.
88        Parameters:
89            none
90        Example:
91            copiedscan = scan.copy()
92        """
93        sd = scantable(sdtable._copy(self))
94        return sd
95
96    def get_scan(self, scanid=None):
97        """
98        Return a specific scan (by scanno) or collection of scans (by
99        source name) in a new scantable.
100        Parameters:
101            scanid:    a scanno or a source name
102        Example:
103            scan.get_scan('323p459')
104            # gets all scans containing the source '323p459'
105        """
106        if scanid is None:
107            print "Please specify a scan no or name to retrieve from the scantable"
108        try:
109            if type(scanid) is str:
110                s = sdtable._getsource(self,scanid)
111                return scantable(s)
112            elif type(scanid) is int:
113                s = sdtable._getscan(self,scanid)
114                return scantable(s)
115        except RuntimeError:
116            print "Couldn't find any match."
117
118    def __str__(self):
119        return sdtable.summary(self)
120
121    def summary(self,filename=None):
122        """
123        Print a summary of the contents of this scantable.
124        Parameters:
125            filename:    the name of a file to write the putput to
126                         Default - no file output
127        """
128        info = sdtable.summary(self)
129        if filename is not None:
130            data = open(filename, 'w')
131            data.write(info)
132            data.close()
133        print info
134
135    def set_selection(self, thebeam=0,theif=0,thepol=0):
136        """
137        Set the spectrum for individual operations.
138        Parameters:
139            thebeam,theif,thepol:    a number
140        Example:
141            scan.set_selection(0,0,1)
142            pol1sig = scan.stats(all=False) # returns std dev for beam=0
143                                            # if=0, pol=1
144        """
145        self.setbeam(thebeam)
146        self.setpol(thepol)
147        self.setif(theif)
148        return
149
150    def get_selection(self):
151        """
152        Return/print a the current 'cursor' into the Beam/IF/Pol cube.
153        Parameters:
154            none
155        Returns:
156            a list of values (currentBeam,currentIF,currentPol)
157        Example:
158            none
159        """
160        i = self.getbeam()
161        j = self.getif()
162        k = self.getpol()
163        if self._vb:
164            print "--------------------------------------------------"
165            print " Cursor selection"
166            print "--------------------------------------------------"
167            out = 'Beam=%d IF=%d Pol=%d '% (i,j,k)
168            print out
169        return i,j,k
170
171    def stats(self, stat='stddev', mask=None, all=True):
172        """
173        Determine the specified statistic of the current beam/if/pol
174        Takes a 'mask' as an optional parameter to specify which
175        channels should be excluded.
176        Parameters:
177            stat:    'min', 'max', 'sumsq', 'sum', 'mean'
178                     'var', 'stddev', 'avdev', 'rms', 'median'
179            mask:    an optional mask specifying where the statistic
180                     should be determined.
181            all:     optional flag to show all (default) or a selected
182                     spectrum of Beam/IF/Pol
183
184        Example:
185            scan.set_unit('channel')
186            msk = scan.create_mask([100,200],[500,600])
187            scan.stats(stat='mean', mask=m)
188        """
189        from asap._asap import stats as _stats
190        if mask == None:
191            mask = ones(self.nchan())
192        if all:
193            out = ''
194            tmp = []
195            for l in range(self.nrow()):
196                tm = self._gettime(l)
197                out += 'Time[%s]:\n' % (tm)
198                for i in range(self.nbeam()):
199                    self.setbeam(i)
200                    if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
201                    for j in range(self.nif()):
202                        self.setif(j)
203                        if self.nif() > 1: out +=  ' IF[%d] ' % (j)
204                        for k in range(self.npol()):
205                            self.setpol(k)
206                            if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
207                            statval = _stats(self,mask,stat)
208                            tmp.append(statval)
209                            out += '= %3.3f\n' % (statval[l])
210                out += "--------------------------------------------------\n"
211
212            if self._vb:
213                print "--------------------------------------------------"
214                print " ",stat
215                out += "--------------------------------------------------\n"
216                print out
217            return tmp
218
219        else:
220            i = self.getbeam()
221            j = self.getif()
222            k = self.getpol()
223            statval = _stats(self,mask,stat)
224            out = ''
225            for l in range(self.nrow()):
226                tm = self._gettime(l)
227                out += 'Time[%s]:\n' % (tm)
228                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
229                if self.nif() > 1: out +=  ' IF[%d] ' % (j)
230                if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
231                out += '= %3.3f\n' % (statval[l])
232                out +=  "--------------------------------------------------\n"
233            if self._vb:
234                print "--------------------------------------------------"
235                print " ",stat
236                print "--------------------------------------------------"
237                print out
238            return statval
239
240    def stddev(self,mask=None, all=True):
241        """
242        Determine the standard deviation of the current beam/if/pol
243        Takes a 'mask' as an optional parameter to specify which
244        channels should be excluded.
245        Parameters:
246            mask:    an optional mask specifying where the standard
247                     deviation should be determined.
248            all:     optional flag to show all or a selected
249                     spectrum of Beam/IF/Pol
250
251        Example:
252            scan.set_unit('channel')
253            msk = scan.create_mask([100,200],[500,600])
254            scan.stddev(mask=m)
255        """
256        return self.stats(stat='stddev',mask=mask, all=all);
257
258    def get_tsys(self, all=True):
259        """
260        Return the System temperatures.
261        Parameters:
262            all:    optional parameter to get the Tsys values for all
263                    Beams/IFs/Pols (default) or just the one selected
264                    with scantable.set_selection()
265                    [True or False]
266        Returns:
267            a list of Tsys values.
268        """
269        if all:
270            out = ''
271            tmp = []
272            for l in range(self.nrow()):
273                tm = self._gettime(l)
274                out += 'Time[%s]:\n' % (tm)
275                for i in range(self.nbeam()):
276                    self.setbeam(i)
277                    if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
278                    for j in range(self.nif()):
279                        self.setif(j)
280                        if self.nif() > 1: out +=  ' IF[%d] ' % (j)
281                        for k in range(self.npol()):
282                            self.setpol(k)
283                            if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
284                            ts = self._gettsys()
285                            tmp.append(ts)
286                            out += '= %3.3f\n' % (ts[l])
287                out+= "--------------------------------------------------\n"
288
289            if self._vb:
290                print "--------------------------------------------------"
291                print " Tsys"
292                print "--------------------------------------------------"
293                print out
294            return tmp
295        else:
296            i = self.getbeam()
297            j = self.getif()
298            k = self.getpol()
299            ts = self._gettsys()
300            out = ''
301            for l in range(self.nrow()):
302                tm = self._gettime(l)
303                out += 'Time[%s]:\n' % (tm)
304                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
305                if self.nif() > 1: out +=  ' IF[%d] ' % (j)
306                if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
307                out += '= %3.3f\n' % (ts[l])
308                out += "--------------------------------------------------\n"
309
310            if self._vb:
311                print "--------------------------------------------------"
312                print " Tsys"
313                print "--------------------------------------------------"
314                print out
315            return ts
316       
317    def get_time(self):
318        """
319        Get a list of time stamps for the observations.
320        Return a string for each integration in the scantable.
321        Parameters:
322            none
323        Example:
324            none
325        """
326        out = []
327        for i in range(self.nrow()):
328            out.append(self._gettime(i))
329        return out
330
331    def set_unit(self, unit='channel'):
332        """
333        Set the unit for all following operations on this scantable
334        Parameters:
335            unit:    optional unit, default is 'channel'
336                     one of '*Hz','km/s','channel', ''
337        """
338
339        if unit in ['','pixel', 'channel']:
340            unit = ''
341        inf = list(self._getcoordinfo())
342        inf[0] = unit
343        self._setcoordinfo(inf)
344        if self._p: self.plot()
345
346    def set_freqframe(self, frame='LSRK'):
347        """
348        Set the frame type of the Spectral Axis.
349        Parameters:
350            frame:   an optional frame type, default 'LSRK'.
351        Examples:
352            scan.set_freqframe('BARY')
353        """
354        valid = ['REST','TOPO','LSRD','LSRK','BARY', \
355                   'GEO','GALACTO','LGROUP','CMB']
356        if frame in valid:
357            inf = list(self._getcoordinfo())
358            inf[1] = frame
359            self._setcoordinfo(inf)
360        else:
361            print "Please specify a valid freq type. Valid types are:\n",valid
362           
363    def get_unit(self):
364        """
365        Get the default unit set in this scantable
366        Parameters:
367        Returns:
368            A unit string
369        """
370        inf = self._getcoordinfo()
371        unit = inf[0]
372        if unit == '': unit = 'channel'
373        return unit
374
375    def get_abcissa(self, rowno=0):
376        """
377        Get the abcissa in the current coordinate setup for the currently
378        selected Beam/IF/Pol
379        Parameters:
380            none
381        Returns:
382            The abcissa values and it's format string.
383        """
384        abc = self.getabcissa(rowno)
385        lbl = self.getabcissalabel(rowno)
386        return abc, lbl
387
388    def create_mask(self, *args, **kwargs):
389        """
390        Compute and return a mask based on [min,max] windows.
391        The specified windows are to be INCLUDED, when the mask is
392        applied.
393        Parameters:
394            [min,max],[min2,max2],...
395                Pairs of start/end points specifying the regions
396                to be masked
397            invert:     optional argument. If specified as True,
398                        return an inverted mask, i.e. the regions
399                        specified are EXCLUDED
400        Example:
401            scan.set_unit('channel')
402
403            a)
404            msk = scan.set_mask([400,500],[800,900])
405            # masks everything outside 400 and 500
406            # and 800 and 900 in the unit 'channel'
407
408            b)
409            msk = scan.set_mask([400,500],[800,900], invert=True)
410            # masks the regions between 400 and 500
411            # and 800 and 900 in the unit 'channel'
412           
413        """
414        u = self._getcoordinfo()[0]
415        if self._vb:
416            if u == "": u = "channel"
417            print "The current mask window unit is", u
418        n = self.nchan()
419        data = self.getabcissa()
420        msk = zeros(n)
421        for  window in args:
422            if (len(window) != 2 or window[0] > window[1] ):
423                print "A window needs to be defined as [min,max]"
424                return
425            for i in range(n):
426                if data[i] >= window[0] and data[i] < window[1]:
427                    msk[i] = 1
428        if kwargs.has_key('invert'):
429            if kwargs.get('invert'):
430                from numarray import logical_not
431                msk = logical_not(msk)
432        return msk
433   
434    def set_restfreqs(self, freqs, unit='Hz'):
435        """
436        Set the restfrequency(s) for this scantable.
437        Parameters:
438            freqs:    one or more frequencies
439            unit:     optional 'unit', default 'Hz'
440        Example:
441            scan.set_restfreqs([1000000000.0])
442        """
443        if type(freqs) is float or int:
444            freqs = (freqs)
445        sdtable._setrestfreqs(self,freqs, unit)
446        return
447
448    def flag_spectrum(self, thebeam, theif, thepol):
449        """
450        This flags a selected spectrum in the scan 'for good'.
451        USE WITH CARE - not reversible.
452        Use masks for non-permanent exclusion of channels.
453        Parameters:
454            thebeam,theif,thepol:    all have to be explicitly
455                                     specified
456        Example:
457            scan.flag_spectrum(0,0,1)
458            flags the spectrum for Beam=0, IF=0, Pol=1
459        """
460        if (thebeam < self.nbeam() and  theif < self.nif() and thepol < self.npol()):
461            stable.setbeam(thebeam)
462            stable.setif(theif)
463            stable.setpol(thepol)
464            stable._flag(self)
465        else:
466            print "Please specify a valid (Beam/IF/Pol)"
467        return
468
469    def plot(self, what='spectrum',col='Pol', panel=None):
470        """
471        Plot the spectra contained in the scan. Alternatively you can also
472        Plot Tsys vs Time
473        Parameters:
474            what:     a choice of 'spectrum' (default) or 'tsys'
475            col:      which out of Beams/IFs/Pols should be colour stacked
476            panel:    set up multiple panels, currently not working.
477        """
478        validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
479
480        validyax = ['spectrum','tsys']
481        from asap.asaplot import ASAPlot
482        if not self._p:
483            self._p = ASAPlot()
484            #print "Plotting not enabled"
485            #return
486        if self._p.is_dead:
487            del self._p
488            self._p = ASAPlot()
489        npan = 1
490        x = None
491        if what == 'tsys':
492            n = self.nrow()
493            if n < 2:
494                print "Only one integration. Can't plot."
495                return
496        self._p.hold()
497        self._p.clear()
498        if panel == 'Time':
499            npan = self.nrow()
500            self._p.set_panels(rows=npan)
501        xlab,ylab,tlab = None,None,None
502        vb = self._verbose()
503        self._verbose(False)
504        sel = self.get_selection()
505        for i in range(npan):
506            self._p.subplot(i)
507            for j in range(validcol[col]):
508                x = None
509                y = None
510                m = None
511                tlab = self._getsourcename(i)
512                import re
513                tlab = re.sub('_S','',tlab)
514                if col == 'Beam':
515                    self.setbeam(j)
516                elif col == 'IF':
517                    self.setif(j)
518                elif col == 'Pol':
519                    self.setpol(j)
520                if what == 'tsys':
521                    x = range(self.nrow())
522                    xlab = 'Time [pixel]'
523                    m = list(ones(len(x)))
524                    y = []
525                    ylab = r'$T_{sys}$'
526                    for k in range(len(x)):
527                        y.append(self._gettsys(k))
528                else:
529                    x,xlab = self.get_abcissa(i)
530                    y = self.getspectrum(i)
531                    ylab = r'Flux'
532                    m = self.getmask(i)
533                llab = col+' '+str(j)
534                self._p.set_line(label=llab)
535                self._p.plot(x,y,m)
536            self._p.set_axes('xlabel',xlab)
537            self._p.set_axes('ylabel',ylab)
538            self._p.set_axes('title',tlab)
539        self._p.release()
540        self.set_selection(sel[0],sel[1],sel[2])
541        self._verbose(vb)
542        return
Note: See TracBrowser for help on using the repository browser.