source: trunk/python/scantable.py@ 343

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

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