source: trunk/python/scantable.py@ 401

Last change on this file since 401 was 393, checked in by kil064, 20 years ago

merge set_restfreqs and select_restfreqs into 1
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.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):
348 """
349 Get a list of time stamps for the observations.
350 Return a string for each integration in the scantable.
351 Parameters:
352 none
353 Example:
354 none
355 """
356 out = []
357 for i in range(self.nrow()):
358 out.append(self._gettime(i))
359 return out
360
361 def set_unit(self, unit='channel'):
362 """
363 Set the unit for all following operations on this scantable
364 Parameters:
365 unit: optional unit, default is 'channel'
366 one of '*Hz','km/s','channel', ''
367 """
368
369 if unit in ['','pixel', 'channel']:
370 unit = ''
371 inf = list(self._getcoordinfo())
372 inf[0] = unit
373 self._setcoordinfo(inf)
374 if self._p: self.plot()
375
376 def set_instrument (self, instr):
377 """
378 Set the instrument for subsequent processing
379 Parameters:
380 instr: Select from "ATPKSMB", "ATPKSHOH", "ATMOPRA",
381 "DSS-43" (Tid), "CEDUNA", and "HOBART"
382 """
383 self._setInstrument(instr)
384
385 def set_doppler(self, doppler='RADIO'):
386 """
387 Set the doppler for all following operations on this scantable.
388 Parameters:
389 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
390 """
391
392 inf = list(self._getcoordinfo())
393 inf[2] = doppler
394 self._setcoordinfo(inf)
395 if self._p: self.plot()
396
397 def set_freqframe(self, frame=None):
398 """
399 Set the frame type of the Spectral Axis.
400 Parameters:
401 frame: an optional frame type, default 'LSRK'.
402 Examples:
403 scan.set_freqframe('BARY')
404 """
405 if not frame: frame = rcParams['scantable.freqframe']
406 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
407 'GEO','GALACTO','LGROUP','CMB']
408 if frame in valid:
409 inf = list(self._getcoordinfo())
410 inf[1] = frame
411 self._setcoordinfo(inf)
412 else:
413 print "Please specify a valid freq type. Valid types are:\n",valid
414
415 def get_unit(self):
416 """
417 Get the default unit set in this scantable
418 Parameters:
419 Returns:
420 A unit string
421 """
422 inf = self._getcoordinfo()
423 unit = inf[0]
424 if unit == '': unit = 'channel'
425 return unit
426
427 def get_abcissa(self, rowno=0):
428 """
429 Get the abcissa in the current coordinate setup for the currently
430 selected Beam/IF/Pol
431 Parameters:
432 rowno: an optional row number in the scantable. Default is the
433 first row, i.e. rowno=0
434 Returns:
435 The abcissa values and it's format string.
436 """
437 abc = self._getabcissa(rowno)
438 lbl = self._getabcissalabel(rowno)
439 return abc, lbl
440
441 def create_mask(self, *args, **kwargs):
442 """
443 Compute and return a mask based on [min,max] windows.
444 The specified windows are to be INCLUDED, when the mask is
445 applied.
446 Parameters:
447 [min,max],[min2,max2],...
448 Pairs of start/end points specifying the regions
449 to be masked
450 invert: optional argument. If specified as True,
451 return an inverted mask, i.e. the regions
452 specified are EXCLUDED
453 Example:
454 scan.set_unit('channel')
455
456 a)
457 msk = scan.set_mask([400,500],[800,900])
458 # masks everything outside 400 and 500
459 # and 800 and 900 in the unit 'channel'
460
461 b)
462 msk = scan.set_mask([400,500],[800,900], invert=True)
463 # masks the regions between 400 and 500
464 # and 800 and 900 in the unit 'channel'
465
466 """
467 u = self._getcoordinfo()[0]
468 if self._vb:
469 if u == "": u = "channel"
470 print "The current mask window unit is", u
471 n = self.nchan()
472 data = self._getabcissa()
473 msk = zeros(n)
474 for window in args:
475 if (len(window) != 2 or window[0] > window[1] ):
476 print "A window needs to be defined as [min,max]"
477 return
478 for i in range(n):
479 if data[i] >= window[0] and data[i] < window[1]:
480 msk[i] = 1
481 if kwargs.has_key('invert'):
482 if kwargs.get('invert'):
483 from numarray import logical_not
484 msk = logical_not(msk)
485 return msk
486
487 def get_restfreqs(self):
488 """
489 Get the restfrequency(s) stored in this scantable.
490 The return value(s) are always of unit 'Hz'
491 Parameters:
492 none
493 Returns:
494 a list of doubles
495 """
496 return list(self._getrestfreqs())
497
498 def set_restfreqs(self, freqs, unit='Hz', source=None, theif=None):
499 """
500 Select the restfrequency for the specified source and IF OR
501 replace for all IFs. If the 'freqs' argument holds a scalar,
502 then that rest frequency will be applied to the selected
503 data (and added to the list of available rest frequencies).
504 In this way, you can set a rest frequency for each
505 source and IF combination. If the 'freqs' argument holds
506 a vector, then it MUST be of length the number of IFs
507 (and the available restfrequencies will be replaced by
508 this vector). In this case, *all* data ('source' and
509 'theif' are ignored) have the restfrequency set per IF according
510 to the corresponding value you give in the 'freqs' vector.
511 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
512 IF 1 gets restfreq 2e9.
513 Parameters:
514 freqs: rest frequencies
515 unit: unit for rest frequency (default 'Hz')
516 source: Source name (blank means all)
517 theif: IF (-1 means all)
518 Example:
519 scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
520 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
521 """
522 if source is None:
523 source = ""
524 if theif is None:
525 theif = -1
526 t = type(freqs)
527 if t is int or t is float:
528 freqs = [freqs]
529 sdtable._setrestfreqs(self, freqs, unit, source, theif)
530 return
531
532
533 def flag_spectrum(self, thebeam, theif, thepol):
534 """
535 This flags a selected spectrum in the scan 'for good'.
536 USE WITH CARE - not reversible.
537 Use masks for non-permanent exclusion of channels.
538 Parameters:
539 thebeam,theif,thepol: all have to be explicitly
540 specified
541 Example:
542 scan.flag_spectrum(0,0,1)
543 flags the spectrum for Beam=0, IF=0, Pol=1
544 """
545 if (thebeam < self.nbeam() and theif < self.nif() and thepol < self.npol()):
546 stable.setbeam(thebeam)
547 stable.setif(theif)
548 stable.setpol(thepol)
549 stable._flag(self)
550 else:
551 print "Please specify a valid (Beam/IF/Pol)"
552 return
553
554 def plot(self, what='spectrum',col='Pol', panel=None):
555 """
556 Plot the spectra contained in the scan. Alternatively you can also
557 Plot Tsys vs Time
558 Parameters:
559 what: a choice of 'spectrum' (default) or 'tsys'
560 col: which out of Beams/IFs/Pols should be colour stacked
561 panel: set up multiple panels, currently not working.
562 """
563 print "Warning! Not fully functional. Use plotter.plot() instead"
564
565 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
566
567 validyax = ['spectrum','tsys']
568 from asap.asaplot import ASAPlot
569 if not self._p:
570 self._p = ASAPlot()
571 #print "Plotting not enabled"
572 #return
573 if self._p.is_dead:
574 del self._p
575 self._p = ASAPlot()
576 npan = 1
577 x = None
578 if what == 'tsys':
579 n = self.nrow()
580 if n < 2:
581 print "Only one integration. Can't plot."
582 return
583 self._p.hold()
584 self._p.clear()
585 if panel == 'Time':
586 npan = self.nrow()
587 self._p.set_panels(rows=npan)
588 xlab,ylab,tlab = None,None,None
589 self._vb = False
590 sel = self.get_cursor()
591 for i in range(npan):
592 if npan > 1:
593 self._p.subplot(i)
594 for j in range(validcol[col]):
595 x = None
596 y = None
597 m = None
598 tlab = self._getsourcename(i)
599 import re
600 tlab = re.sub('_S','',tlab)
601 if col == 'Beam':
602 self.setbeam(j)
603 elif col == 'IF':
604 self.setif(j)
605 elif col == 'Pol':
606 self.setpol(j)
607 if what == 'tsys':
608 x = range(self.nrow())
609 xlab = 'Time [pixel]'
610 m = list(ones(len(x)))
611 y = []
612 ylab = r'$T_{sys}$'
613 for k in range(len(x)):
614 y.append(self._gettsys(k))
615 else:
616 x,xlab = self.get_abcissa(i)
617 y = self._getspectrum(i)
618 ylab = r'Flux'
619 m = self._getmask(i)
620 llab = col+' '+str(j)
621 self._p.set_line(label=llab)
622 self._p.plot(x,y,m)
623 self._p.set_axes('xlabel',xlab)
624 self._p.set_axes('ylabel',ylab)
625 self._p.set_axes('title',tlab)
626 self._p.release()
627 self.set_cursor(sel[0],sel[1],sel[2])
628 self._vb = rcParams['verbose']
629 return
630
631 print out
632
633 def _print_values(self, dat, label='', timestamps=[]):
634 d = dat['data']
635 a = dat['axes']
636 shp = d.getshape()
637 out = ''
638 for i in range(shp[3]):
639 out += '%s [%s]:\n' % (a[3],timestamps[i])
640 t = d[:,:,:,i]
641 for j in range(shp[0]):
642 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
643 for k in range(shp[1]):
644 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
645 for l in range(shp[2]):
646 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
647 out += '= %3.3f\n' % (t[j,k,l])
648 out += "--------------------------------------------------\n"
649 print "--------------------------------------------------"
650 print " ", label
651 print "--------------------------------------------------"
652 print out
Note: See TracBrowser for help on using the repository browser.