source: trunk/python/scantable.py@ 442

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

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