source: trunk/python/scantable.py@ 218

Last change on this file since 218 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
RevLine 
[102]1from asap._asap import sdtable
[189]2from numarray import ones,zeros
[102]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
[196]51 def save(self, name="", format='ASAP'):
[116]52 """
53 Store the scantable on disk. This can be a asap file or SDFITS/MS2.
54 Parameters:
[196]55 name: the name of the outputfile. For format="FITS" this
56 is the root file name (and can be empty).
[116]57 format: an optional file format. Default is ASAP.
[194]58 Allowed are - 'ASAP' (save as ASAP Table),
59 'SDFITS' (save as SDFITS file)
60 'FITS' (saves each row as a FITS Image)
[200]61 'ASCII' (saves as ascii text file)
[194]62 'MS2' (saves as an aips++ MeasurementSet V2)
[116]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
[194]71 w = _sw(format)
72 w.write(self, name)
[116]73 return
74
[113]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 """
[126]80 if len(args) == 0:
81 return self._vb
82 elif type(args[0]) is bool:
[113]83 self._vb = args[0]
84 return
85
[102]86 def copy(self):
87 """
88 Return a copy of this scantable.
89 Parameters:
[113]90 none
[102]91 Example:
92 copiedscan = scan.copy()
93 """
[113]94 sd = scantable(sdtable._copy(self))
95 return sd
96
[102]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:
[113]104 scan.get_scan('323p459')
[102]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:
[113]111 s = sdtable._getsource(self,scanid)
[102]112 return scantable(s)
113 elif type(scanid) is int:
[113]114 s = sdtable._getscan(self,scanid)
[102]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)
[135]143 pol1sig = scan.stats(all=False) # returns std dev for beam=0
[181]144 # if=0, pol=1
[102]145 """
146 self.setbeam(thebeam)
147 self.setpol(thepol)
148 self.setif(theif)
149 return
150
151 def get_selection(self):
[113]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 """
[102]161 i = self.getbeam()
162 j = self.getif()
163 k = self.getpol()
[113]164 if self._vb:
[181]165 print "--------------------------------------------------"
166 print " Cursor selection"
167 print "--------------------------------------------------"
[113]168 out = 'Beam=%d IF=%d Pol=%d '% (i,j,k)
169 print out
[102]170 return i,j,k
171
[135]172 def stats(self, stat='stddev', mask=None, all=True):
[102]173 """
[135]174 Determine the specified statistic of the current beam/if/pol
[102]175 Takes a 'mask' as an optional parameter to specify which
176 channels should be excluded.
177 Parameters:
[135]178 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
179 'var', 'stddev', 'avdev', 'rms', 'median'
180 mask: an optional mask specifying where the statistic
[102]181 should be determined.
[156]182 all: optional flag to show all (default) or a selected
[102]183 spectrum of Beam/IF/Pol
184
185 Example:
[113]186 scan.set_unit('channel')
[102]187 msk = scan.create_mask([100,200],[500,600])
[135]188 scan.stats(stat='mean', mask=m)
[102]189 """
[135]190 from asap._asap import stats as _stats
[102]191 if mask == None:
192 mask = ones(self.nchan())
193 if all:
194 out = ''
195 tmp = []
[181]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
[102]218 return tmp
219
220 else:
221 i = self.getbeam()
222 j = self.getif()
223 k = self.getpol()
[181]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
[102]240
[135]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
[102]259 def get_tsys(self, all=True):
[113]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 """
[102]270 if all:
[181]271 out = ''
[102]272 tmp = []
[181]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
[102]295 return tmp
296 else:
297 i = self.getbeam()
298 j = self.getif()
299 k = self.getpol()
[113]300 ts = self._gettsys()
[181]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
[102]316 return ts
[113]317
318 def get_time(self):
319 """
320 Get a list of time stamps for the observations.
[181]321 Return a string for each integration in the scantable.
[113]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
[102]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'
[113]337 one of '*Hz','km/s','channel', ''
[102]338 """
[113]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)
[102]361 else:
[113]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
[102]375
[158]376 def get_abcissa(self, rowno=0):
[102]377 """
[158]378 Get the abcissa in the current coordinate setup for the currently
[113]379 selected Beam/IF/Pol
380 Parameters:
381 none
382 Returns:
[158]383 The abcissa values and it's format string.
[113]384 """
[158]385 abc = self.getabcissa(rowno)
386 lbl = self.getabcissalabel(rowno)
387 return abc, lbl
[113]388
389 def create_mask(self, *args, **kwargs):
390 """
[102]391 Compute and return a mask based on [min,max] windows.
[189]392 The specified windows are to be INCLUDED, when the mask is
[113]393 applied.
[102]394 Parameters:
395 [min,max],[min2,max2],...
396 Pairs of start/end points specifying the regions
397 to be masked
[189]398 invert: optional argument. If specified as True,
399 return an inverted mask, i.e. the regions
400 specified are EXCLUDED
[102]401 Example:
[113]402 scan.set_unit('channel')
403
404 a)
[102]405 msk = scan.set_mask([400,500],[800,900])
[189]406 # masks everything outside 400 and 500
[113]407 # and 800 and 900 in the unit 'channel'
408
409 b)
410 msk = scan.set_mask([400,500],[800,900], invert=True)
[189]411 # masks the regions between 400 and 500
[113]412 # and 800 and 900 in the unit 'channel'
413
[102]414 """
[113]415 u = self._getcoordinfo()[0]
416 if self._vb:
417 if u == "": u = "channel"
418 print "The current mask window unit is", u
[102]419 n = self.nchan()
[158]420 data = self.getabcissa()
[189]421 msk = zeros(n)
[102]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):
[113]427 if data[i] >= window[0] and data[i] < window[1]:
[189]428 msk[i] = 1
[113]429 if kwargs.has_key('invert'):
430 if kwargs.get('invert'):
431 from numarray import logical_not
432 msk = logical_not(msk)
[102]433 return msk
[113]434
435 def set_restfreqs(self, freqs, unit='Hz'):
[102]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:
[113]442 scan.set_restfreqs([1000000000.0])
[102]443 """
[113]444 if type(freqs) is float or int:
445 freqs = (freqs)
446 sdtable._setrestfreqs(self,freqs, unit)
[102]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
[113]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'
[135]476 col: which out of Beams/IFs/Pols should be colour stacked
[113]477 panel: set up multiple panels, currently not working.
478 """
479 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
[124]480
[113]481 validyax = ['spectrum','tsys']
[189]482 from asap.asaplot import ASAPlot
[113]483 if not self._p:
484 self._p = ASAPlot()
[189]485 #print "Plotting not enabled"
486 #return
487 if self._p.is_dead:
488 del self._p
489 self._p = ASAPlot()
[113]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
[124]497 self._p.hold()
498 self._p.clear()
[113]499 if panel == 'Time':
500 npan = self.nrow()
[124]501 self._p.set_panels(rows=npan)
[113]502 xlab,ylab,tlab = None,None,None
[124]503 vb = self._verbose()
[113]504 self._verbose(False)
505 sel = self.get_selection()
506 for i in range(npan):
[124]507 self._p.subplot(i)
[113]508 for j in range(validcol[col]):
509 x = None
510 y = None
511 m = None
512 tlab = self._getsourcename(i)
[124]513 import re
514 tlab = re.sub('_S','',tlab)
[113]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:
[158]530 x,xlab = self.get_abcissa(i)
[113]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)
[124]537 self._p.set_axes('xlabel',xlab)
538 self._p.set_axes('ylabel',ylab)
539 self._p.set_axes('title',tlab)
[113]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.