source: trunk/python/scantable.py@ 418

Last change on this file since 418 was 415, checked in by mar637, 20 years ago

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