source: trunk/python/scantable.py@ 304

Last change on this file since 304 was 283, checked in by mar637, 20 years ago

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