source: trunk/python/scantable.py@ 459

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

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