source: trunk/python/scantable.py@ 476

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

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