source: trunk/python/scantable.py@ 182

Last change on this file since 182 was 181, checked in by mar637, 20 years ago

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