source: trunk/python/scantable.py@ 402

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

add lines arg to funciton set_restfreqs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.5 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:
[381]141 s = sdtable._getscan(self,[scanid])
142 return scantable(s)
143 elif type(scanid) is list:
[113]144 s = sdtable._getscan(self,scanid)
[102]145 return scantable(s)
[381]146 else:
147 print "Illegal scanid type, use 'int' or 'list' if ints."
[102]148 except RuntimeError:
149 print "Couldn't find any match."
150
151 def __str__(self):
[381]152 return sdtable._summary(self,True)
[102]153
[381]154 def summary(self,filename=None, verbose=None):
[102]155 """
156 Print a summary of the contents of this scantable.
157 Parameters:
158 filename: the name of a file to write the putput to
159 Default - no file output
[381]160 verbose: print extra info such as the frequency table
161 The default (False) is taken from .asaprc
[102]162 """
[381]163 info = sdtable._summary(self, verbose)
164 if verbose is None: verbose = rcParams['scantable.verbosesummary']
[102]165 if filename is not None:
[256]166 if filename is "":
167 filename = 'scantable_summary.txt'
[102]168 data = open(filename, 'w')
169 data.write(info)
170 data.close()
171 print info
172
[256]173 def set_cursor(self, thebeam=0,theif=0,thepol=0):
[102]174 """
175 Set the spectrum for individual operations.
176 Parameters:
177 thebeam,theif,thepol: a number
178 Example:
[256]179 scan.set_cursor(0,0,1)
[135]180 pol1sig = scan.stats(all=False) # returns std dev for beam=0
[181]181 # if=0, pol=1
[102]182 """
183 self.setbeam(thebeam)
184 self.setpol(thepol)
185 self.setif(theif)
186 return
187
[256]188 def get_cursor(self):
[113]189 """
190 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
191 Parameters:
192 none
193 Returns:
194 a list of values (currentBeam,currentIF,currentPol)
195 Example:
196 none
197 """
[102]198 i = self.getbeam()
199 j = self.getif()
200 k = self.getpol()
[113]201 if self._vb:
[181]202 print "--------------------------------------------------"
[256]203 print " Cursor position"
[181]204 print "--------------------------------------------------"
[226]205 out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
[113]206 print out
[102]207 return i,j,k
208
[226]209 def stats(self, stat='stddev', mask=None, all=None):
[102]210 """
[135]211 Determine the specified statistic of the current beam/if/pol
[102]212 Takes a 'mask' as an optional parameter to specify which
213 channels should be excluded.
214 Parameters:
[135]215 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
216 'var', 'stddev', 'avdev', 'rms', 'median'
217 mask: an optional mask specifying where the statistic
[102]218 should be determined.
[241]219 all: if true show all (default or .asaprc) rather
220 that the cursor selected spectrum of Beam/IF/Pol
[102]221
222 Example:
[113]223 scan.set_unit('channel')
[102]224 msk = scan.create_mask([100,200],[500,600])
[135]225 scan.stats(stat='mean', mask=m)
[102]226 """
[226]227 if all is None: all = rcParams['scantable.allaxes']
[135]228 from asap._asap import stats as _stats
[256]229 from numarray import array,zeros,Float
[102]230 if mask == None:
231 mask = ones(self.nchan())
[256]232 axes = ['Beam','IF','Pol','Time']
233
[102]234 if all:
[256]235 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
236 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
237 arr = array(zeros(n),shape=shp,type=Float)
238
239 for i in range(self.nbeam()):
240 self.setbeam(i)
241 for j in range(self.nif()):
242 self.setif(j)
243 for k in range(self.npol()):
244 self.setpol(k)
245 arr[i,j,k,:] = _stats(self,mask,stat,-1)
246 retval = {'axes': axes, 'data': arr, 'cursor':None}
247 tm = [self._gettime(val) for val in range(self.nrow())]
[181]248 if self._vb:
[256]249 self._print_values(retval,stat,tm)
250 return retval
[102]251
252 else:
[256]253 i,j,k = (self.getbeam(),self.getif(),self.getpol())
[241]254 statval = _stats(self,mask,stat,-1)
[181]255 out = ''
256 for l in range(self.nrow()):
257 tm = self._gettime(l)
258 out += 'Time[%s]:\n' % (tm)
259 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
260 if self.nif() > 1: out += ' IF[%d] ' % (j)
261 if self.npol() > 1: out += ' Pol[%d] ' % (k)
262 out += '= %3.3f\n' % (statval[l])
263 out += "--------------------------------------------------\n"
[226]264
[181]265 if self._vb:
266 print "--------------------------------------------------"
267 print " ",stat
268 print "--------------------------------------------------"
269 print out
[256]270 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
271 return retval
[102]272
[226]273 def stddev(self,mask=None, all=None):
[135]274 """
275 Determine the standard deviation of the current beam/if/pol
276 Takes a 'mask' as an optional parameter to specify which
277 channels should be excluded.
278 Parameters:
279 mask: an optional mask specifying where the standard
280 deviation should be determined.
[226]281 all: optional flag to show all or a cursor selected
282 spectrum of Beam/IF/Pol. Default is all or taken
283 from .asaprc
[135]284
285 Example:
286 scan.set_unit('channel')
287 msk = scan.create_mask([100,200],[500,600])
288 scan.stddev(mask=m)
289 """
[226]290 if all is None: all = rcParams['scantable.allaxes']
[135]291 return self.stats(stat='stddev',mask=mask, all=all);
292
[226]293 def get_tsys(self, all=None):
[113]294 """
295 Return the System temperatures.
296 Parameters:
297 all: optional parameter to get the Tsys values for all
298 Beams/IFs/Pols (default) or just the one selected
[256]299 with scantable.set_cursor()
[113]300 [True or False]
301 Returns:
302 a list of Tsys values.
303 """
[226]304 if all is None: all = rcParams['scantable.allaxes']
[256]305 from numarray import array,zeros,Float
306 axes = ['Beam','IF','Pol','Time']
307
[102]308 if all:
[256]309 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
310 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
311 arr = array(zeros(n),shape=shp,type=Float)
312
313 for i in range(self.nbeam()):
314 self.setbeam(i)
315 for j in range(self.nif()):
316 self.setif(j)
317 for k in range(self.npol()):
318 self.setpol(k)
319 arr[i,j,k,:] = self._gettsys()
320 retval = {'axes': axes, 'data': arr, 'cursor':None}
321 tm = [self._gettime(val) for val in range(self.nrow())]
[181]322 if self._vb:
[256]323 self._print_values(retval,'Tsys',tm)
324 return retval
325
[102]326 else:
[256]327 i,j,k = (self.getbeam(),self.getif(),self.getpol())
328 statval = self._gettsys()
[181]329 out = ''
330 for l in range(self.nrow()):
331 tm = self._gettime(l)
332 out += 'Time[%s]:\n' % (tm)
333 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
334 if self.nif() > 1: out += ' IF[%d] ' % (j)
335 if self.npol() > 1: out += ' Pol[%d] ' % (k)
[256]336 out += '= %3.3f\n' % (statval[l])
337 out += "--------------------------------------------------\n"
338
[181]339 if self._vb:
340 print "--------------------------------------------------"
[256]341 print " TSys"
[181]342 print "--------------------------------------------------"
343 print out
[256]344 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
345 return retval
[113]346
347 def get_time(self):
348 """
349 Get a list of time stamps for the observations.
[181]350 Return a string for each integration in the scantable.
[113]351 Parameters:
352 none
353 Example:
354 none
355 """
356 out = []
357 for i in range(self.nrow()):
358 out.append(self._gettime(i))
359 return out
[102]360
361 def set_unit(self, unit='channel'):
362 """
363 Set the unit for all following operations on this scantable
364 Parameters:
365 unit: optional unit, default is 'channel'
[113]366 one of '*Hz','km/s','channel', ''
[102]367 """
[113]368
369 if unit in ['','pixel', 'channel']:
370 unit = ''
371 inf = list(self._getcoordinfo())
372 inf[0] = unit
373 self._setcoordinfo(inf)
374 if self._p: self.plot()
375
[358]376 def set_instrument (self, instr):
377 """
378 Set the instrument for subsequent processing
379 Parameters:
380 instr: Select from "ATPKSMB", "ATPKSHOH", "ATMOPRA",
381 "DSS-43" (Tid), "CEDUNA", and "HOBART"
382 """
383 self._setInstrument(instr)
384
[276]385 def set_doppler(self, doppler='RADIO'):
386 """
387 Set the doppler for all following operations on this scantable.
388 Parameters:
389 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
390 """
391
392 inf = list(self._getcoordinfo())
393 inf[2] = doppler
394 self._setcoordinfo(inf)
395 if self._p: self.plot()
396
[226]397 def set_freqframe(self, frame=None):
[113]398 """
399 Set the frame type of the Spectral Axis.
400 Parameters:
401 frame: an optional frame type, default 'LSRK'.
402 Examples:
403 scan.set_freqframe('BARY')
404 """
[226]405 if not frame: frame = rcParams['scantable.freqframe']
[113]406 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
407 'GEO','GALACTO','LGROUP','CMB']
408 if frame in valid:
409 inf = list(self._getcoordinfo())
410 inf[1] = frame
411 self._setcoordinfo(inf)
[102]412 else:
[113]413 print "Please specify a valid freq type. Valid types are:\n",valid
414
415 def get_unit(self):
416 """
417 Get the default unit set in this scantable
418 Parameters:
419 Returns:
420 A unit string
421 """
422 inf = self._getcoordinfo()
423 unit = inf[0]
424 if unit == '': unit = 'channel'
425 return unit
[102]426
[158]427 def get_abcissa(self, rowno=0):
[102]428 """
[158]429 Get the abcissa in the current coordinate setup for the currently
[113]430 selected Beam/IF/Pol
431 Parameters:
[226]432 rowno: an optional row number in the scantable. Default is the
433 first row, i.e. rowno=0
[113]434 Returns:
[158]435 The abcissa values and it's format string.
[113]436 """
[256]437 abc = self._getabcissa(rowno)
438 lbl = self._getabcissalabel(rowno)
[158]439 return abc, lbl
[113]440
441 def create_mask(self, *args, **kwargs):
442 """
[102]443 Compute and return a mask based on [min,max] windows.
[189]444 The specified windows are to be INCLUDED, when the mask is
[113]445 applied.
[102]446 Parameters:
447 [min,max],[min2,max2],...
448 Pairs of start/end points specifying the regions
449 to be masked
[189]450 invert: optional argument. If specified as True,
451 return an inverted mask, i.e. the regions
452 specified are EXCLUDED
[102]453 Example:
[113]454 scan.set_unit('channel')
455
456 a)
[102]457 msk = scan.set_mask([400,500],[800,900])
[189]458 # masks everything outside 400 and 500
[113]459 # and 800 and 900 in the unit 'channel'
460
461 b)
462 msk = scan.set_mask([400,500],[800,900], invert=True)
[189]463 # masks the regions between 400 and 500
[113]464 # and 800 and 900 in the unit 'channel'
465
[102]466 """
[113]467 u = self._getcoordinfo()[0]
468 if self._vb:
469 if u == "": u = "channel"
470 print "The current mask window unit is", u
[102]471 n = self.nchan()
[256]472 data = self._getabcissa()
[189]473 msk = zeros(n)
[102]474 for window in args:
475 if (len(window) != 2 or window[0] > window[1] ):
476 print "A window needs to be defined as [min,max]"
477 return
478 for i in range(n):
[113]479 if data[i] >= window[0] and data[i] < window[1]:
[189]480 msk[i] = 1
[113]481 if kwargs.has_key('invert'):
482 if kwargs.get('invert'):
483 from numarray import logical_not
484 msk = logical_not(msk)
[102]485 return msk
[113]486
[256]487 def get_restfreqs(self):
488 """
489 Get the restfrequency(s) stored in this scantable.
490 The return value(s) are always of unit 'Hz'
491 Parameters:
492 none
493 Returns:
494 a list of doubles
495 """
496 return list(self._getrestfreqs())
[102]497
[402]498 def lines(self):
[391]499 """
[402]500 Print the list of known spectral lines
501 """
502 sdtable._lines(self)
503
504 def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None, theif=None):
505 """
[393]506 Select the restfrequency for the specified source and IF OR
507 replace for all IFs. If the 'freqs' argument holds a scalar,
508 then that rest frequency will be applied to the selected
509 data (and added to the list of available rest frequencies).
510 In this way, you can set a rest frequency for each
511 source and IF combination. If the 'freqs' argument holds
512 a vector, then it MUST be of length the number of IFs
513 (and the available restfrequencies will be replaced by
514 this vector). In this case, *all* data ('source' and
515 'theif' are ignored) have the restfrequency set per IF according
516 to the corresponding value you give in the 'freqs' vector.
517 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
518 IF 1 gets restfreq 2e9.
[402]519
520 You can also specify the frequencies via known line names
521 in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
522 takes precedence. See the list of known names via function
523 scantable.lines()
[391]524 Parameters:
[402]525 freqs: list of rest frequencies
[393]526 unit: unit for rest frequency (default 'Hz')
[402]527 lines: list of known spectral lines names (alternative to freqs).
528 See possible list via scantable.lines()
[391]529 source: Source name (blank means all)
530 theif: IF (-1 means all)
531 Example:
[393]532 scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
533 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
[391]534 """
535 if source is None:
536 source = ""
537 if theif is None:
538 theif = -1
[393]539 t = type(freqs)
540 if t is int or t is float:
541 freqs = [freqs]
[402]542 if freqs is None:
543 freqs = []
544 t = type(lines)
545 if t is str:
546 lines = [lines]
547 if lines is None:
548 lines = []
549 sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
[391]550 return
551
552
[102]553 def flag_spectrum(self, thebeam, theif, thepol):
554 """
555 This flags a selected spectrum in the scan 'for good'.
556 USE WITH CARE - not reversible.
557 Use masks for non-permanent exclusion of channels.
558 Parameters:
559 thebeam,theif,thepol: all have to be explicitly
560 specified
561 Example:
562 scan.flag_spectrum(0,0,1)
563 flags the spectrum for Beam=0, IF=0, Pol=1
564 """
565 if (thebeam < self.nbeam() and theif < self.nif() and thepol < self.npol()):
566 stable.setbeam(thebeam)
567 stable.setif(theif)
568 stable.setpol(thepol)
569 stable._flag(self)
570 else:
571 print "Please specify a valid (Beam/IF/Pol)"
572 return
[113]573
574 def plot(self, what='spectrum',col='Pol', panel=None):
575 """
576 Plot the spectra contained in the scan. Alternatively you can also
577 Plot Tsys vs Time
578 Parameters:
579 what: a choice of 'spectrum' (default) or 'tsys'
[135]580 col: which out of Beams/IFs/Pols should be colour stacked
[113]581 panel: set up multiple panels, currently not working.
582 """
[283]583 print "Warning! Not fully functional. Use plotter.plot() instead"
584
[113]585 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
[124]586
[113]587 validyax = ['spectrum','tsys']
[189]588 from asap.asaplot import ASAPlot
[113]589 if not self._p:
590 self._p = ASAPlot()
[189]591 #print "Plotting not enabled"
592 #return
593 if self._p.is_dead:
594 del self._p
595 self._p = ASAPlot()
[113]596 npan = 1
597 x = None
598 if what == 'tsys':
599 n = self.nrow()
600 if n < 2:
601 print "Only one integration. Can't plot."
602 return
[124]603 self._p.hold()
604 self._p.clear()
[113]605 if panel == 'Time':
606 npan = self.nrow()
[124]607 self._p.set_panels(rows=npan)
[113]608 xlab,ylab,tlab = None,None,None
[226]609 self._vb = False
[256]610 sel = self.get_cursor()
[113]611 for i in range(npan):
[226]612 if npan > 1:
613 self._p.subplot(i)
[113]614 for j in range(validcol[col]):
615 x = None
616 y = None
617 m = None
618 tlab = self._getsourcename(i)
[124]619 import re
620 tlab = re.sub('_S','',tlab)
[113]621 if col == 'Beam':
622 self.setbeam(j)
623 elif col == 'IF':
624 self.setif(j)
625 elif col == 'Pol':
626 self.setpol(j)
627 if what == 'tsys':
628 x = range(self.nrow())
629 xlab = 'Time [pixel]'
630 m = list(ones(len(x)))
631 y = []
632 ylab = r'$T_{sys}$'
633 for k in range(len(x)):
634 y.append(self._gettsys(k))
635 else:
[158]636 x,xlab = self.get_abcissa(i)
[256]637 y = self._getspectrum(i)
[113]638 ylab = r'Flux'
[256]639 m = self._getmask(i)
[113]640 llab = col+' '+str(j)
641 self._p.set_line(label=llab)
642 self._p.plot(x,y,m)
[124]643 self._p.set_axes('xlabel',xlab)
644 self._p.set_axes('ylabel',ylab)
645 self._p.set_axes('title',tlab)
[113]646 self._p.release()
[256]647 self.set_cursor(sel[0],sel[1],sel[2])
[226]648 self._vb = rcParams['verbose']
[113]649 return
[256]650
651 print out
652
653 def _print_values(self, dat, label='', timestamps=[]):
654 d = dat['data']
655 a = dat['axes']
656 shp = d.getshape()
657 out = ''
658 for i in range(shp[3]):
659 out += '%s [%s]:\n' % (a[3],timestamps[i])
660 t = d[:,:,:,i]
661 for j in range(shp[0]):
662 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
663 for k in range(shp[1]):
664 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
665 for l in range(shp[2]):
666 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
667 out += '= %3.3f\n' % (t[j,k,l])
668 out += "--------------------------------------------------\n"
669 print "--------------------------------------------------"
670 print " ", label
671 print "--------------------------------------------------"
672 print out
Note: See TracBrowser for help on using the repository browser.