source: trunk/python/scantable.py @ 195

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

allow FITS format in function 'save'

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