source: trunk/python/scantable.py @ 415

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

Made note about reader limitation on rpfits.
fixed up summary (again)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.1 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, unit=None):
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           unit:         brightness unit; must be consistent with K or Jy.
24                         Over-rides the default selected by the reader
25                         (input rpfits/sdfits/ms) or replaces the value
26                         in existing scantables
27        """
28        self._vb = rcParams['verbose']
29        self._p = None
30        from os import stat as st
31        import stat
32        if isinstance(filename,sdtable):
33            sdtable.__init__(self, filename)           
34            if unit is not None:
35                self.set_fluxunit(unit)                       
36        else:
37            import os.path
38            if not os.path.exists(filename):
39                print "File '%s' not found." % (filename)
40                return
41            filename = os.path.expandvars(filename)
42            if os.path.isdir(filename):
43                # crude check if asap table
44                if os.path.exists(filename+'/table.info'):
45                    sdtable.__init__(self, filename)
46                    if unit is not None:
47                        self.set_fluxunit(unit)                       
48                else:
49                    print "The given file '%s'is not a valid asap table." % (filename)
50                    return
51            else:
52                autoav = rcParams['scantable.autoaverage']
53
54                from asap._asap import sdreader
55                ifSel = -1
56                beamSel = -1
57                r = sdreader(filename,ifSel,beamSel)
58                print 'Importing data...'
59                r._read([-1])
60                tbl = r._getdata()
61                if unit is not None:
62                    tbl.set_fluxunit(unit)
63                if autoav:
64                    from asap._asap import average
65                    tmp = tuple([tbl])
66                    print 'Auto averaging integrations...'
67                    tbl2 = average(tmp,(),True,'none')
68                    sdtable.__init__(self,tbl2)
69                    del r,tbl
70                else:
71                    sdtable.__init__(self,tbl)
72
73    def save(self, name=None, format=None, overwrite=False):
74        """
75        Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
76        Image FITS or MS2 format.
77        Parameters:
78            name:        the name of the outputfile. For format="FITS" this
79                         is the directory file name into which all the files will
80                         be written (default is 'asap_FITS')
81            format:      an optional file format. Default is ASAP.
82                         Allowed are - 'ASAP' (save as ASAP [aips++] Table),
83                                       'SDFITS' (save as SDFITS file)
84                                       'FITS' (saves each row as a FITS Image)
85                                       'ASCII' (saves as ascii text file)
86                                       'MS2' (saves as an aips++
87                                              MeasurementSet V2)
88            overwrite:   If the file should be overwritten if it exists.
89                         The default False is to return with warning
90                         without writing the output. USE WITH CARE.
91        Example:
92            scan.save('myscan.asap')
93            scan.save('myscan.sdfits','SDFITS')
94        """
95        from os import path
96        if format is None: format = rcParams['scantable.save']
97        suffix = '.'+format.lower()
98        if name is None or name =="":
99            name = 'scantable'+suffix
100            print "No filename given. Using default name %s..." % name
101        name = path.expandvars(name)
102        if path.isfile(name) or path.isdir(name):
103            if not overwrite:
104                print "File %s already exists." % name
105                return
106        if format == 'ASAP':
107            self._save(name)
108        else:
109            from asap._asap import sdwriter as _sw
110            w = _sw(format)
111            w.write(self, name)
112        return
113
114    def copy(self):
115        """
116        Return a copy of this scantable.
117        Parameters:
118            none
119        Example:
120            copiedscan = scan.copy()
121        """
122        sd = scantable(sdtable._copy(self))
123        return sd
124
125    def get_scan(self, scanid=None):
126        """
127        Return a specific scan (by scanno) or collection of scans (by
128        source name) in a new scantable.
129        Parameters:
130            scanid:    a scanno or a source name
131        Example:
132            scan.get_scan('323p459')
133            # gets all scans containing the source '323p459'
134        """
135        if scanid is None:
136            print "Please specify a scan no or name to retrieve from the scantable"
137        try:
138            if type(scanid) is str:
139                s = sdtable._getsource(self,scanid)
140                return scantable(s)
141            elif type(scanid) is int:
142                s = sdtable._getscan(self,[scanid])
143                return scantable(s)
144            elif type(scanid) is list:
145                s = sdtable._getscan(self,scanid)
146                return scantable(s)
147            else:
148                print "Illegal scanid type, use 'int' or 'list' if ints."
149        except RuntimeError:
150            print "Couldn't find any match."
151
152    def __str__(self):
153        return sdtable._summary(self,True)
154
155    def summary(self,filename=None, verbose=None):
156        """
157        Print a summary of the contents of this scantable.
158        Parameters:
159            filename:    the name of a file to write the putput to
160                         Default - no file output
161            verbose:     print extra info such as the frequency table
162                         The default (False) is taken from .asaprc
163        """
164        info = sdtable._summary(self, verbose)
165        if verbose is None: verbose = rcParams['scantable.verbosesummary']
166        if filename is not None:
167            if filename is "":
168                filename = 'scantable_summary.txt'
169            from os.path import expandvars, isdir
170            filename = expandvars(filename)
171            if not isdir(filename):
172                data = open(filename, 'w')
173                data.write(info)
174                data.close()
175            else:
176                print "Illegal file name '%s'." % (filename)
177        print info
178
179    def set_cursor(self, thebeam=0,theif=0,thepol=0):
180        """
181        Set the spectrum for individual operations.
182        Parameters:
183            thebeam,theif,thepol:    a number
184        Example:
185            scan.set_cursor(0,0,1)
186            pol1sig = scan.stats(all=False) # returns std dev for beam=0
187                                            # if=0, pol=1
188        """
189        self.setbeam(thebeam)
190        self.setpol(thepol)
191        self.setif(theif)
192        return
193
194    def get_cursor(self):
195        """
196        Return/print a the current 'cursor' into the Beam/IF/Pol cube.
197        Parameters:
198            none
199        Returns:
200            a list of values (currentBeam,currentIF,currentPol)
201        Example:
202            none
203        """
204        i = self.getbeam()
205        j = self.getif()
206        k = self.getpol()
207        if self._vb:
208            print "--------------------------------------------------"
209            print " Cursor position"
210            print "--------------------------------------------------"
211            out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
212            print out
213        return i,j,k
214
215    def stats(self, stat='stddev', mask=None, all=None):
216        """
217        Determine the specified statistic of the current beam/if/pol
218        Takes a 'mask' as an optional parameter to specify which
219        channels should be excluded.
220        Parameters:
221            stat:    'min', 'max', 'sumsq', 'sum', 'mean'
222                     'var', 'stddev', 'avdev', 'rms', 'median'
223            mask:    an optional mask specifying where the statistic
224                     should be determined.
225            all:     if true show all (default or .asaprc) rather
226                     that the cursor selected spectrum of Beam/IF/Pol
227
228        Example:
229            scan.set_unit('channel')
230            msk = scan.create_mask([100,200],[500,600])
231            scan.stats(stat='mean', mask=m)
232        """
233        if all is None: all = rcParams['scantable.allaxes']
234        from asap._asap import stats as _stats
235        from numarray import array,zeros,Float
236        if mask == None:
237            mask = ones(self.nchan())
238        axes = ['Beam','IF','Pol','Time']
239
240        if all:
241            n = self.nbeam()*self.nif()*self.npol()*self.nrow()
242            shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
243            arr = array(zeros(n),shape=shp,type=Float)
244
245            for i in range(self.nbeam()):
246                self.setbeam(i)
247                for j in range(self.nif()):
248                    self.setif(j)
249                    for k in range(self.npol()):
250                        self.setpol(k)
251                        arr[i,j,k,:] = _stats(self,mask,stat,-1)
252            retval = {'axes': axes, 'data': arr, 'cursor':None}
253            tm = [self._gettime(val) for val in range(self.nrow())]
254            if self._vb:
255                self._print_values(retval,stat,tm)
256            return retval
257
258        else:
259            i,j,k = (self.getbeam(),self.getif(),self.getpol())
260            statval = _stats(self,mask,stat,-1)
261            out = ''
262            for l in range(self.nrow()):
263                tm = self._gettime(l)
264                out += 'Time[%s]:\n' % (tm)
265                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
266                if self.nif() > 1: out +=  ' IF[%d] ' % (j)
267                if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
268                out += '= %3.3f\n' % (statval[l])
269                out +=  "--------------------------------------------------\n"
270
271            if self._vb:
272                print "--------------------------------------------------"
273                print " ",stat
274                print "--------------------------------------------------"
275                print out
276            retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
277            return retval
278
279    def stddev(self,mask=None, all=None):
280        """
281        Determine the standard deviation of the current beam/if/pol
282        Takes a 'mask' as an optional parameter to specify which
283        channels should be excluded.
284        Parameters:
285            mask:    an optional mask specifying where the standard
286                     deviation should be determined.
287            all:     optional flag to show all or a cursor selected
288                     spectrum of Beam/IF/Pol. Default is all or taken
289                     from .asaprc
290
291        Example:
292            scan.set_unit('channel')
293            msk = scan.create_mask([100,200],[500,600])
294            scan.stddev(mask=m)
295        """
296        if all is None: all = rcParams['scantable.allaxes']
297        return self.stats(stat='stddev',mask=mask, all=all);
298
299    def get_tsys(self, all=None):
300        """
301        Return the System temperatures.
302        Parameters:
303            all:    optional parameter to get the Tsys values for all
304                    Beams/IFs/Pols (default) or just the one selected
305                    with scantable.set_cursor()
306                    [True or False]
307        Returns:
308            a list of Tsys values.
309        """
310        if all is None: all = rcParams['scantable.allaxes']
311        from numarray import array,zeros,Float
312        axes = ['Beam','IF','Pol','Time']
313
314        if all:
315            n = self.nbeam()*self.nif()*self.npol()*self.nrow()
316            shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
317            arr = array(zeros(n),shape=shp,type=Float)
318
319            for i in range(self.nbeam()):
320                self.setbeam(i)
321                for j in range(self.nif()):
322                    self.setif(j)
323                    for k in range(self.npol()):
324                        self.setpol(k)
325                        arr[i,j,k,:] = self._gettsys()
326            retval = {'axes': axes, 'data': arr, 'cursor':None}
327            tm = [self._gettime(val) for val in range(self.nrow())]
328            if self._vb:
329                self._print_values(retval,'Tsys',tm)
330            return retval
331
332        else:
333            i,j,k = (self.getbeam(),self.getif(),self.getpol())
334            statval = self._gettsys()
335            out = ''
336            for l in range(self.nrow()):
337                tm = self._gettime(l)
338                out += 'Time[%s]:\n' % (tm)
339                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
340                if self.nif() > 1: out +=  ' IF[%d] ' % (j)
341                if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
342                out += '= %3.3f\n' % (statval[l])
343                out +=  "--------------------------------------------------\n"
344
345            if self._vb:
346                print "--------------------------------------------------"
347                print " TSys"
348                print "--------------------------------------------------"
349                print out
350            retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
351            return retval
352       
353    def get_time(self, row=-1):
354        """
355        Get a list of time stamps for the observations.
356        Return a string for each integration in the scantable.
357        Parameters:
358            row:    row no of integration. Default -1 return all rows
359        Example:
360            none
361        """
362        out = []
363        if row == -1:
364            for i in range(self.nrow()):
365                out.append(self._gettime(i))
366            return out
367        else:
368            if row < self.nrow():
369                return self._gettime(row)
370
371    def set_unit(self, unit='channel'):
372        """
373        Set the unit for all following operations on this scantable
374        Parameters:
375            unit:    optional unit, default is 'channel'
376                     one of '*Hz','km/s','channel', ''
377        """
378
379        if unit in ['','pixel', 'channel']:
380            unit = ''
381        inf = list(self._getcoordinfo())
382        inf[0] = unit
383        self._setcoordinfo(inf)
384        if self._p: self.plot()
385
386    def set_instrument (self, instr):
387        """
388        Set the instrument for subsequent processing
389        Parameters:
390            instr:    Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
391                      'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
392        """
393        self._setInstrument(instr)
394
395    def set_doppler(self, doppler='RADIO'):
396        """
397        Set the doppler for all following operations on this scantable.
398        Parameters:
399            doppler:    One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
400        """
401
402        inf = list(self._getcoordinfo())
403        inf[2] = doppler
404        self._setcoordinfo(inf)
405        if self._p: self.plot()
406
407    def set_freqframe(self, frame=None):
408        """
409        Set the frame type of the Spectral Axis.
410        Parameters:
411            frame:   an optional frame type, default 'LSRK'.
412        Examples:
413            scan.set_freqframe('BARY')
414        """
415        if not frame: frame = rcParams['scantable.freqframe']
416        valid = ['REST','TOPO','LSRD','LSRK','BARY', \
417                   'GEO','GALACTO','LGROUP','CMB']
418        if frame in valid:
419            inf = list(self._getcoordinfo())
420            inf[1] = frame
421            self._setcoordinfo(inf)
422        else:
423            print "Please specify a valid freq type. Valid types are:\n",valid
424           
425    def get_unit(self):
426        """
427        Get the default unit set in this scantable
428        Parameters:
429        Returns:
430            A unit string
431        """
432        inf = self._getcoordinfo()
433        unit = inf[0]
434        if unit == '': unit = 'channel'
435        return unit
436
437    def get_abcissa(self, rowno=0):
438        """
439        Get the abcissa in the current coordinate setup for the currently
440        selected Beam/IF/Pol
441        Parameters:
442            rowno:    an optional row number in the scantable. Default is the
443                      first row, i.e. rowno=0
444        Returns:
445            The abcissa values and it's format string (as a dictionary)
446        """
447        abc = self._getabcissa(rowno)
448        lbl = self._getabcissalabel(rowno)       
449        return abc, lbl
450        #return {'abcissa':abc,'label':lbl}
451
452    def create_mask(self, *args, **kwargs):
453        """
454        Compute and return a mask based on [min,max] windows.
455        The specified windows are to be INCLUDED, when the mask is
456        applied.
457        Parameters:
458            [min,max],[min2,max2],...
459                Pairs of start/end points specifying the regions
460                to be masked
461            invert:     optional argument. If specified as True,
462                        return an inverted mask, i.e. the regions
463                        specified are EXCLUDED
464        Example:
465            scan.set_unit('channel')
466
467            a)
468            msk = scan.set_mask([400,500],[800,900])
469            # masks everything outside 400 and 500
470            # and 800 and 900 in the unit 'channel'
471
472            b)
473            msk = scan.set_mask([400,500],[800,900], invert=True)
474            # masks the regions between 400 and 500
475            # and 800 and 900 in the unit 'channel'
476           
477        """
478        u = self._getcoordinfo()[0]
479        if self._vb:
480            if u == "": u = "channel"
481            print "The current mask window unit is", u
482        n = self.nchan()
483        data = self._getabcissa()
484        msk = zeros(n)
485        for  window in args:
486            if (len(window) != 2 or window[0] > window[1] ):
487                print "A window needs to be defined as [min,max]"
488                return
489            for i in range(n):
490                if data[i] >= window[0] and data[i] < window[1]:
491                    msk[i] = 1
492        if kwargs.has_key('invert'):
493            if kwargs.get('invert'):
494                from numarray import logical_not
495                msk = logical_not(msk)
496        return msk
497   
498    def get_restfreqs(self):
499        """
500        Get the restfrequency(s) stored in this scantable.
501        The return value(s) are always of unit 'Hz'
502        Parameters:
503            none
504        Returns:
505            a list of doubles
506        """
507        return list(self._getrestfreqs())
508
509    def lines(self):
510        """
511        Print the list of known spectral lines
512        """
513        sdtable._lines(self)
514
515    def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None, theif=None):
516        """
517        Select the restfrequency for the specified source and IF OR
518        replace for all IFs.  If the 'freqs' argument holds a scalar,
519        then that rest frequency will be applied to the selected
520        data (and added to the list of available rest frequencies).
521        In this way, you can set a rest frequency for each
522        source and IF combination.   If the 'freqs' argument holds
523        a vector, then it MUST be of length the number of IFs
524        (and the available restfrequencies will be replaced by
525        this vector).  In this case, *all* data ('source' and
526        'theif' are ignored) have the restfrequency set per IF according
527        to the corresponding value you give in the 'freqs' vector. 
528        E.g. 'freqs=[1e9,2e9]'  would mean IF 0 gets restfreq 1e9 and
529        IF 1 gets restfreq 2e9.
530
531        You can also specify the frequencies via known line names
532        in the argument 'lines'.  Use 'freqs' or 'lines'.  'freqs'
533        takes precedence. See the list of known names via function
534        scantable.lines()
535        Parameters:
536            freqs:   list of rest frequencies
537            unit:    unit for rest frequency (default 'Hz')
538            lines:   list of known spectral lines names (alternative to freqs).
539                     See possible list via scantable.lines()
540            source:  Source name (blank means all)
541            theif:   IF (-1 means all)
542        Example:
543            scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
544            scan.set_restfreqs(freqs=[1.4e9,1.67e9])
545        """
546        if source is None:
547            source = ""
548        if theif is None:
549            theif = -1
550        t = type(freqs)
551        if t is int or t is float:
552           freqs = [freqs]
553        if freqs is None:
554           freqs = []
555        t = type(lines)
556        if t is str:
557           lines = [lines]
558        if lines is None:
559           lines = []
560        sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
561        return
562
563
564    def flag_spectrum(self, thebeam, theif, thepol):
565        """
566        This flags a selected spectrum in the scan 'for good'.
567        USE WITH CARE - not reversible.
568        Use masks for non-permanent exclusion of channels.
569        Parameters:
570            thebeam,theif,thepol:    all have to be explicitly
571                                     specified
572        Example:
573            scan.flag_spectrum(0,0,1)
574            flags the spectrum for Beam=0, IF=0, Pol=1
575        """
576        if (thebeam < self.nbeam() and
577            theif < self.nif() and
578            thepol < self.npol()):
579            sdtable.setbeam(self, thebeam)
580            sdtable.setif(self, theif)
581            sdtable.setpol(self, thepol)
582            sdtable._flag(self)
583        else:
584            print "Please specify a valid (Beam/IF/Pol)"
585        return
586
587    def plot(self, what='spectrum',col='Pol', panel=None):
588        """
589        Plot the spectra contained in the scan. Alternatively you can also
590        Plot Tsys vs Time
591        Parameters:
592            what:     a choice of 'spectrum' (default) or 'tsys'
593            col:      which out of Beams/IFs/Pols should be colour stacked
594            panel:    set up multiple panels, currently not working.
595        """
596        print "Warning! Not fully functional. Use plotter.plot() instead"
597       
598        validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
599
600        validyax = ['spectrum','tsys']
601        from asap.asaplot import ASAPlot
602        if not self._p:
603            self._p = ASAPlot()
604            #print "Plotting not enabled"
605            #return
606        if self._p.is_dead:
607            del self._p
608            self._p = ASAPlot()
609        npan = 1
610        x = None
611        if what == 'tsys':
612            n = self.nrow()
613            if n < 2:
614                print "Only one integration. Can't plot."
615                return
616        self._p.hold()
617        self._p.clear()
618        if panel == 'Time':
619            npan = self.nrow()
620            self._p.set_panels(rows=npan)
621        xlab,ylab,tlab = None,None,None
622        self._vb = False
623        sel = self.get_cursor()       
624        for i in range(npan):
625            if npan > 1:
626                self._p.subplot(i)
627            for j in range(validcol[col]):
628                x = None
629                y = None
630                m = None
631                tlab = self._getsourcename(i)
632                import re
633                tlab = re.sub('_S','',tlab)
634                if col == 'Beam':
635                    self.setbeam(j)
636                elif col == 'IF':
637                    self.setif(j)
638                elif col == 'Pol':
639                    self.setpol(j)
640                if what == 'tsys':
641                    x = range(self.nrow())
642                    xlab = 'Time [pixel]'
643                    m = list(ones(len(x)))
644                    y = []
645                    ylab = r'$T_{sys}$'
646                    for k in range(len(x)):
647                        y.append(self._gettsys(k))
648                else:
649                    x,xlab = self.get_abcissa(i)
650                    y = self._getspectrum(i)
651                    ylab = r'Flux'
652                    m = self._getmask(i)
653                llab = col+' '+str(j)
654                self._p.set_line(label=llab)
655                self._p.plot(x,y,m)
656            self._p.set_axes('xlabel',xlab)
657            self._p.set_axes('ylabel',ylab)
658            self._p.set_axes('title',tlab)
659        self._p.release()
660        self.set_cursor(sel[0],sel[1],sel[2])
661        self._vb = rcParams['verbose']
662        return
663
664        print out
665
666    def _print_values(self, dat, label='', timestamps=[]):
667        d = dat['data']
668        a = dat['axes']
669        shp = d.getshape()
670        out = ''
671        for i in range(shp[3]):
672            out += '%s [%s]:\n' % (a[3],timestamps[i])
673            t = d[:,:,:,i]
674            for j in range(shp[0]):
675                if shp[0] > 1: out +=  ' %s[%d] ' % (a[0],j)
676                for k in range(shp[1]):
677                    if shp[1] > 1: out +=  ' %s[%d] ' % (a[1],k)
678                    for l in range(shp[2]):
679                        if shp[2] > 1: out +=  ' %s[%d] ' % (a[2],l)
680                        out += '= %3.3f\n' % (t[j,k,l])
681            out += "--------------------------------------------------\n"
682        print "--------------------------------------------------"
683        print " ", label
684        print "--------------------------------------------------"
685        print out
Note: See TracBrowser for help on using the repository browser.