source: trunk/python/scantable.py@ 220

Last change on this file since 220 was 200, checked in by kil064, 20 years ago

document ascii file save capability

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