source: trunk/python/scantable.py@ 144

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

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