source: trunk/python/scantable.py@ 359

Last change on this file since 359 was 358, checked in by kil064, 20 years ago

add binding to function set_instrument (instead of direct pythin binding)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.3 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:
[346]35 self.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
[345]57 r = sdreader(filename,ifSel,beamSel)
[181]58 print 'Importing data...'
59 r.read([-1])
60 tbl = r.getdata()
[346]61 if unit is not None:
62 tbl.set_fluxunit(unit)
[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
[358]368 def set_instrument (self, instr):
369 """
370 Set the instrument for subsequent processing
371 Parameters:
372 instr: Select from "ATPKSMB", "ATPKSHOH", "ATMOPRA",
373 "DSS-43" (Tid), "CEDUNA", and "HOBART"
374 """
375 self._setInstrument(instr)
376
[276]377 def set_doppler(self, doppler='RADIO'):
378 """
379 Set the doppler for all following operations on this scantable.
380 Parameters:
381 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
382 """
383
384 inf = list(self._getcoordinfo())
385 inf[2] = doppler
386 self._setcoordinfo(inf)
387 if self._p: self.plot()
388
[226]389 def set_freqframe(self, frame=None):
[113]390 """
391 Set the frame type of the Spectral Axis.
392 Parameters:
393 frame: an optional frame type, default 'LSRK'.
394 Examples:
395 scan.set_freqframe('BARY')
396 """
[226]397 if not frame: frame = rcParams['scantable.freqframe']
[113]398 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
399 'GEO','GALACTO','LGROUP','CMB']
400 if frame in valid:
401 inf = list(self._getcoordinfo())
402 inf[1] = frame
403 self._setcoordinfo(inf)
[102]404 else:
[113]405 print "Please specify a valid freq type. Valid types are:\n",valid
406
407 def get_unit(self):
408 """
409 Get the default unit set in this scantable
410 Parameters:
411 Returns:
412 A unit string
413 """
414 inf = self._getcoordinfo()
415 unit = inf[0]
416 if unit == '': unit = 'channel'
417 return unit
[102]418
[158]419 def get_abcissa(self, rowno=0):
[102]420 """
[158]421 Get the abcissa in the current coordinate setup for the currently
[113]422 selected Beam/IF/Pol
423 Parameters:
[226]424 rowno: an optional row number in the scantable. Default is the
425 first row, i.e. rowno=0
[113]426 Returns:
[158]427 The abcissa values and it's format string.
[113]428 """
[256]429 abc = self._getabcissa(rowno)
430 lbl = self._getabcissalabel(rowno)
[158]431 return abc, lbl
[113]432
433 def create_mask(self, *args, **kwargs):
434 """
[102]435 Compute and return a mask based on [min,max] windows.
[189]436 The specified windows are to be INCLUDED, when the mask is
[113]437 applied.
[102]438 Parameters:
439 [min,max],[min2,max2],...
440 Pairs of start/end points specifying the regions
441 to be masked
[189]442 invert: optional argument. If specified as True,
443 return an inverted mask, i.e. the regions
444 specified are EXCLUDED
[102]445 Example:
[113]446 scan.set_unit('channel')
447
448 a)
[102]449 msk = scan.set_mask([400,500],[800,900])
[189]450 # masks everything outside 400 and 500
[113]451 # and 800 and 900 in the unit 'channel'
452
453 b)
454 msk = scan.set_mask([400,500],[800,900], invert=True)
[189]455 # masks the regions between 400 and 500
[113]456 # and 800 and 900 in the unit 'channel'
457
[102]458 """
[113]459 u = self._getcoordinfo()[0]
460 if self._vb:
461 if u == "": u = "channel"
462 print "The current mask window unit is", u
[102]463 n = self.nchan()
[256]464 data = self._getabcissa()
[189]465 msk = zeros(n)
[102]466 for window in args:
467 if (len(window) != 2 or window[0] > window[1] ):
468 print "A window needs to be defined as [min,max]"
469 return
470 for i in range(n):
[113]471 if data[i] >= window[0] and data[i] < window[1]:
[189]472 msk[i] = 1
[113]473 if kwargs.has_key('invert'):
474 if kwargs.get('invert'):
475 from numarray import logical_not
476 msk = logical_not(msk)
[102]477 return msk
[113]478
479 def set_restfreqs(self, freqs, unit='Hz'):
[102]480 """
481 Set the restfrequency(s) for this scantable.
482 Parameters:
483 freqs: one or more frequencies
484 unit: optional 'unit', default 'Hz'
485 Example:
[113]486 scan.set_restfreqs([1000000000.0])
[102]487 """
[113]488 if type(freqs) is float or int:
489 freqs = (freqs)
490 sdtable._setrestfreqs(self,freqs, unit)
[102]491 return
[256]492 def get_restfreqs(self):
493 """
494 Get the restfrequency(s) stored in this scantable.
495 The return value(s) are always of unit 'Hz'
496 Parameters:
497 none
498 Returns:
499 a list of doubles
500 """
501 return list(self._getrestfreqs())
[102]502
503 def flag_spectrum(self, thebeam, theif, thepol):
504 """
505 This flags a selected spectrum in the scan 'for good'.
506 USE WITH CARE - not reversible.
507 Use masks for non-permanent exclusion of channels.
508 Parameters:
509 thebeam,theif,thepol: all have to be explicitly
510 specified
511 Example:
512 scan.flag_spectrum(0,0,1)
513 flags the spectrum for Beam=0, IF=0, Pol=1
514 """
515 if (thebeam < self.nbeam() and theif < self.nif() and thepol < self.npol()):
516 stable.setbeam(thebeam)
517 stable.setif(theif)
518 stable.setpol(thepol)
519 stable._flag(self)
520 else:
521 print "Please specify a valid (Beam/IF/Pol)"
522 return
[113]523
524 def plot(self, what='spectrum',col='Pol', panel=None):
525 """
526 Plot the spectra contained in the scan. Alternatively you can also
527 Plot Tsys vs Time
528 Parameters:
529 what: a choice of 'spectrum' (default) or 'tsys'
[135]530 col: which out of Beams/IFs/Pols should be colour stacked
[113]531 panel: set up multiple panels, currently not working.
532 """
[283]533 print "Warning! Not fully functional. Use plotter.plot() instead"
534
[113]535 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
[124]536
[113]537 validyax = ['spectrum','tsys']
[189]538 from asap.asaplot import ASAPlot
[113]539 if not self._p:
540 self._p = ASAPlot()
[189]541 #print "Plotting not enabled"
542 #return
543 if self._p.is_dead:
544 del self._p
545 self._p = ASAPlot()
[113]546 npan = 1
547 x = None
548 if what == 'tsys':
549 n = self.nrow()
550 if n < 2:
551 print "Only one integration. Can't plot."
552 return
[124]553 self._p.hold()
554 self._p.clear()
[113]555 if panel == 'Time':
556 npan = self.nrow()
[124]557 self._p.set_panels(rows=npan)
[113]558 xlab,ylab,tlab = None,None,None
[226]559 self._vb = False
[256]560 sel = self.get_cursor()
[113]561 for i in range(npan):
[226]562 if npan > 1:
563 self._p.subplot(i)
[113]564 for j in range(validcol[col]):
565 x = None
566 y = None
567 m = None
568 tlab = self._getsourcename(i)
[124]569 import re
570 tlab = re.sub('_S','',tlab)
[113]571 if col == 'Beam':
572 self.setbeam(j)
573 elif col == 'IF':
574 self.setif(j)
575 elif col == 'Pol':
576 self.setpol(j)
577 if what == 'tsys':
578 x = range(self.nrow())
579 xlab = 'Time [pixel]'
580 m = list(ones(len(x)))
581 y = []
582 ylab = r'$T_{sys}$'
583 for k in range(len(x)):
584 y.append(self._gettsys(k))
585 else:
[158]586 x,xlab = self.get_abcissa(i)
[256]587 y = self._getspectrum(i)
[113]588 ylab = r'Flux'
[256]589 m = self._getmask(i)
[113]590 llab = col+' '+str(j)
591 self._p.set_line(label=llab)
592 self._p.plot(x,y,m)
[124]593 self._p.set_axes('xlabel',xlab)
594 self._p.set_axes('ylabel',ylab)
595 self._p.set_axes('title',tlab)
[113]596 self._p.release()
[256]597 self.set_cursor(sel[0],sel[1],sel[2])
[226]598 self._vb = rcParams['verbose']
[113]599 return
[256]600
601 print out
602
603 def _print_values(self, dat, label='', timestamps=[]):
604 d = dat['data']
605 a = dat['axes']
606 shp = d.getshape()
607 out = ''
608 for i in range(shp[3]):
609 out += '%s [%s]:\n' % (a[3],timestamps[i])
610 t = d[:,:,:,i]
611 for j in range(shp[0]):
612 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
613 for k in range(shp[1]):
614 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
615 for l in range(shp[2]):
616 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
617 out += '= %3.3f\n' % (t[j,k,l])
618 out += "--------------------------------------------------\n"
619 print "--------------------------------------------------"
620 print " ", label
621 print "--------------------------------------------------"
622 print out
Note: See TracBrowser for help on using the repository browser.