source: trunk/python/scantable.py@ 241

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

modify stats function to pass row number arg to C++ statistics

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