source: trunk/python/scantable.py@ 340

Last change on this file since 340 was 340, checked in by kil064, 21 years ago

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