source: trunk/python/scantable.py@ 395

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

merge set_restfreqs and select_restfreqs into 1
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.8 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
[393]498 def set_restfreqs(self, freqs, unit='Hz', source=None, theif=None):
[391]499 """
[393]500 Select the restfrequency for the specified source and IF OR
501 replace for all IFs. If the 'freqs' argument holds a scalar,
502 then that rest frequency will be applied to the selected
503 data (and added to the list of available rest frequencies).
504 In this way, you can set a rest frequency for each
505 source and IF combination. If the 'freqs' argument holds
506 a vector, then it MUST be of length the number of IFs
507 (and the available restfrequencies will be replaced by
508 this vector). In this case, *all* data ('source' and
509 'theif' are ignored) have the restfrequency set per IF according
510 to the corresponding value you give in the 'freqs' vector.
511 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
512 IF 1 gets restfreq 2e9.
[391]513 Parameters:
[393]514 freqs: rest frequencies
515 unit: unit for rest frequency (default 'Hz')
[391]516 source: Source name (blank means all)
517 theif: IF (-1 means all)
518 Example:
[393]519 scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
520 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
[391]521 """
522 if source is None:
523 source = ""
524 if theif is None:
525 theif = -1
[393]526 t = type(freqs)
527 if t is int or t is float:
528 freqs = [freqs]
529 sdtable._setrestfreqs(self, freqs, unit, source, theif)
[391]530 return
531
532
[102]533 def flag_spectrum(self, thebeam, theif, thepol):
534 """
535 This flags a selected spectrum in the scan 'for good'.
536 USE WITH CARE - not reversible.
537 Use masks for non-permanent exclusion of channels.
538 Parameters:
539 thebeam,theif,thepol: all have to be explicitly
540 specified
541 Example:
542 scan.flag_spectrum(0,0,1)
543 flags the spectrum for Beam=0, IF=0, Pol=1
544 """
545 if (thebeam < self.nbeam() and theif < self.nif() and thepol < self.npol()):
546 stable.setbeam(thebeam)
547 stable.setif(theif)
548 stable.setpol(thepol)
549 stable._flag(self)
550 else:
551 print "Please specify a valid (Beam/IF/Pol)"
552 return
[113]553
554 def plot(self, what='spectrum',col='Pol', panel=None):
555 """
556 Plot the spectra contained in the scan. Alternatively you can also
557 Plot Tsys vs Time
558 Parameters:
559 what: a choice of 'spectrum' (default) or 'tsys'
[135]560 col: which out of Beams/IFs/Pols should be colour stacked
[113]561 panel: set up multiple panels, currently not working.
562 """
[283]563 print "Warning! Not fully functional. Use plotter.plot() instead"
564
[113]565 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
[124]566
[113]567 validyax = ['spectrum','tsys']
[189]568 from asap.asaplot import ASAPlot
[113]569 if not self._p:
570 self._p = ASAPlot()
[189]571 #print "Plotting not enabled"
572 #return
573 if self._p.is_dead:
574 del self._p
575 self._p = ASAPlot()
[113]576 npan = 1
577 x = None
578 if what == 'tsys':
579 n = self.nrow()
580 if n < 2:
581 print "Only one integration. Can't plot."
582 return
[124]583 self._p.hold()
584 self._p.clear()
[113]585 if panel == 'Time':
586 npan = self.nrow()
[124]587 self._p.set_panels(rows=npan)
[113]588 xlab,ylab,tlab = None,None,None
[226]589 self._vb = False
[256]590 sel = self.get_cursor()
[113]591 for i in range(npan):
[226]592 if npan > 1:
593 self._p.subplot(i)
[113]594 for j in range(validcol[col]):
595 x = None
596 y = None
597 m = None
598 tlab = self._getsourcename(i)
[124]599 import re
600 tlab = re.sub('_S','',tlab)
[113]601 if col == 'Beam':
602 self.setbeam(j)
603 elif col == 'IF':
604 self.setif(j)
605 elif col == 'Pol':
606 self.setpol(j)
607 if what == 'tsys':
608 x = range(self.nrow())
609 xlab = 'Time [pixel]'
610 m = list(ones(len(x)))
611 y = []
612 ylab = r'$T_{sys}$'
613 for k in range(len(x)):
614 y.append(self._gettsys(k))
615 else:
[158]616 x,xlab = self.get_abcissa(i)
[256]617 y = self._getspectrum(i)
[113]618 ylab = r'Flux'
[256]619 m = self._getmask(i)
[113]620 llab = col+' '+str(j)
621 self._p.set_line(label=llab)
622 self._p.plot(x,y,m)
[124]623 self._p.set_axes('xlabel',xlab)
624 self._p.set_axes('ylabel',ylab)
625 self._p.set_axes('title',tlab)
[113]626 self._p.release()
[256]627 self.set_cursor(sel[0],sel[1],sel[2])
[226]628 self._vb = rcParams['verbose']
[113]629 return
[256]630
631 print out
632
633 def _print_values(self, dat, label='', timestamps=[]):
634 d = dat['data']
635 a = dat['axes']
636 shp = d.getshape()
637 out = ''
638 for i in range(shp[3]):
639 out += '%s [%s]:\n' % (a[3],timestamps[i])
640 t = d[:,:,:,i]
641 for j in range(shp[0]):
642 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
643 for k in range(shp[1]):
644 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
645 for l in range(shp[2]):
646 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
647 out += '= %3.3f\n' % (t[j,k,l])
648 out += "--------------------------------------------------\n"
649 print "--------------------------------------------------"
650 print " ", label
651 print "--------------------------------------------------"
652 print out
Note: See TracBrowser for help on using the repository browser.