source: trunk/python/scantable.py@ 282

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

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