source: trunk/python/scantable.py@ 194

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

allow FITS format in function 'save'

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