source: trunk/python/scantable.py@ 388

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