source: trunk/python/scantable.py @ 200

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

document ascii file save capability

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