source: trunk/python/scantable.py @ 226

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

added rcParams to support rc style default parameters, read from .asaprc

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