source: trunk/python/scantable.py@ 408

Last change on this file since 408 was 407, checked in by mar637, 20 years ago

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