source: trunk/python/scantable.py@ 198

Last change on this file since 198 was 196, checked in by kil064, 20 years ago

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