source: trunk/python/scantable.py @ 189

Last change on this file since 189 was 189, checked in by mar637, 19 years ago

* empty log message *

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