source: trunk/python/scantable.py@ 277

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

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