source: trunk/python/scantable.py@ 373

Last change on this file since 373 was 358, checked in by kil064, 20 years ago

add binding to function set_instrument (instead of direct pythin binding)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.3 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 except RuntimeError:
144 print "Couldn't find any match."
145
146 def __str__(self):
147 return sdtable._summary(self)
148
149 def summary(self,filename=None):
150 """
151 Print a summary of the contents of this scantable.
152 Parameters:
153 filename: the name of a file to write the putput to
154 Default - no file output
155 """
156 info = sdtable._summary(self)
157 if filename is not None:
158 if filename is "":
159 filename = 'scantable_summary.txt'
160 data = open(filename, 'w')
161 data.write(info)
162 data.close()
163 print info
164
165 def set_cursor(self, thebeam=0,theif=0,thepol=0):
166 """
167 Set the spectrum for individual operations.
168 Parameters:
169 thebeam,theif,thepol: a number
170 Example:
171 scan.set_cursor(0,0,1)
172 pol1sig = scan.stats(all=False) # returns std dev for beam=0
173 # if=0, pol=1
174 """
175 self.setbeam(thebeam)
176 self.setpol(thepol)
177 self.setif(theif)
178 return
179
180 def get_cursor(self):
181 """
182 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
183 Parameters:
184 none
185 Returns:
186 a list of values (currentBeam,currentIF,currentPol)
187 Example:
188 none
189 """
190 i = self.getbeam()
191 j = self.getif()
192 k = self.getpol()
193 if self._vb:
194 print "--------------------------------------------------"
195 print " Cursor position"
196 print "--------------------------------------------------"
197 out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
198 print out
199 return i,j,k
200
201 def stats(self, stat='stddev', mask=None, all=None):
202 """
203 Determine the specified statistic of the current beam/if/pol
204 Takes a 'mask' as an optional parameter to specify which
205 channels should be excluded.
206 Parameters:
207 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
208 'var', 'stddev', 'avdev', 'rms', 'median'
209 mask: an optional mask specifying where the statistic
210 should be determined.
211 all: if true show all (default or .asaprc) rather
212 that the cursor selected spectrum of Beam/IF/Pol
213
214 Example:
215 scan.set_unit('channel')
216 msk = scan.create_mask([100,200],[500,600])
217 scan.stats(stat='mean', mask=m)
218 """
219 if all is None: all = rcParams['scantable.allaxes']
220 from asap._asap import stats as _stats
221 from numarray import array,zeros,Float
222 if mask == None:
223 mask = ones(self.nchan())
224 axes = ['Beam','IF','Pol','Time']
225
226 if all:
227 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
228 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
229 arr = array(zeros(n),shape=shp,type=Float)
230
231 for i in range(self.nbeam()):
232 self.setbeam(i)
233 for j in range(self.nif()):
234 self.setif(j)
235 for k in range(self.npol()):
236 self.setpol(k)
237 arr[i,j,k,:] = _stats(self,mask,stat,-1)
238 retval = {'axes': axes, 'data': arr, 'cursor':None}
239 tm = [self._gettime(val) for val in range(self.nrow())]
240 if self._vb:
241 self._print_values(retval,stat,tm)
242 return retval
243
244 else:
245 i,j,k = (self.getbeam(),self.getif(),self.getpol())
246 statval = _stats(self,mask,stat,-1)
247 out = ''
248 for l in range(self.nrow()):
249 tm = self._gettime(l)
250 out += 'Time[%s]:\n' % (tm)
251 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
252 if self.nif() > 1: out += ' IF[%d] ' % (j)
253 if self.npol() > 1: out += ' Pol[%d] ' % (k)
254 out += '= %3.3f\n' % (statval[l])
255 out += "--------------------------------------------------\n"
256
257 if self._vb:
258 print "--------------------------------------------------"
259 print " ",stat
260 print "--------------------------------------------------"
261 print out
262 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
263 return retval
264
265 def stddev(self,mask=None, all=None):
266 """
267 Determine the standard deviation of the current beam/if/pol
268 Takes a 'mask' as an optional parameter to specify which
269 channels should be excluded.
270 Parameters:
271 mask: an optional mask specifying where the standard
272 deviation should be determined.
273 all: optional flag to show all or a cursor selected
274 spectrum of Beam/IF/Pol. Default is all or taken
275 from .asaprc
276
277 Example:
278 scan.set_unit('channel')
279 msk = scan.create_mask([100,200],[500,600])
280 scan.stddev(mask=m)
281 """
282 if all is None: all = rcParams['scantable.allaxes']
283 return self.stats(stat='stddev',mask=mask, all=all);
284
285 def get_tsys(self, all=None):
286 """
287 Return the System temperatures.
288 Parameters:
289 all: optional parameter to get the Tsys values for all
290 Beams/IFs/Pols (default) or just the one selected
291 with scantable.set_cursor()
292 [True or False]
293 Returns:
294 a list of Tsys values.
295 """
296 if all is None: all = rcParams['scantable.allaxes']
297 from numarray import array,zeros,Float
298 axes = ['Beam','IF','Pol','Time']
299
300 if all:
301 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
302 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
303 arr = array(zeros(n),shape=shp,type=Float)
304
305 for i in range(self.nbeam()):
306 self.setbeam(i)
307 for j in range(self.nif()):
308 self.setif(j)
309 for k in range(self.npol()):
310 self.setpol(k)
311 arr[i,j,k,:] = self._gettsys()
312 retval = {'axes': axes, 'data': arr, 'cursor':None}
313 tm = [self._gettime(val) for val in range(self.nrow())]
314 if self._vb:
315 self._print_values(retval,'Tsys',tm)
316 return retval
317
318 else:
319 i,j,k = (self.getbeam(),self.getif(),self.getpol())
320 statval = self._gettsys()
321 out = ''
322 for l in range(self.nrow()):
323 tm = self._gettime(l)
324 out += 'Time[%s]:\n' % (tm)
325 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
326 if self.nif() > 1: out += ' IF[%d] ' % (j)
327 if self.npol() > 1: out += ' Pol[%d] ' % (k)
328 out += '= %3.3f\n' % (statval[l])
329 out += "--------------------------------------------------\n"
330
331 if self._vb:
332 print "--------------------------------------------------"
333 print " TSys"
334 print "--------------------------------------------------"
335 print out
336 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
337 return retval
338
339 def get_time(self):
340 """
341 Get a list of time stamps for the observations.
342 Return a string for each integration in the scantable.
343 Parameters:
344 none
345 Example:
346 none
347 """
348 out = []
349 for i in range(self.nrow()):
350 out.append(self._gettime(i))
351 return out
352
353 def set_unit(self, unit='channel'):
354 """
355 Set the unit for all following operations on this scantable
356 Parameters:
357 unit: optional unit, default is 'channel'
358 one of '*Hz','km/s','channel', ''
359 """
360
361 if unit in ['','pixel', 'channel']:
362 unit = ''
363 inf = list(self._getcoordinfo())
364 inf[0] = unit
365 self._setcoordinfo(inf)
366 if self._p: self.plot()
367
368 def set_instrument (self, instr):
369 """
370 Set the instrument for subsequent processing
371 Parameters:
372 instr: Select from "ATPKSMB", "ATPKSHOH", "ATMOPRA",
373 "DSS-43" (Tid), "CEDUNA", and "HOBART"
374 """
375 self._setInstrument(instr)
376
377 def set_doppler(self, doppler='RADIO'):
378 """
379 Set the doppler for all following operations on this scantable.
380 Parameters:
381 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
382 """
383
384 inf = list(self._getcoordinfo())
385 inf[2] = doppler
386 self._setcoordinfo(inf)
387 if self._p: self.plot()
388
389 def set_freqframe(self, frame=None):
390 """
391 Set the frame type of the Spectral Axis.
392 Parameters:
393 frame: an optional frame type, default 'LSRK'.
394 Examples:
395 scan.set_freqframe('BARY')
396 """
397 if not frame: frame = rcParams['scantable.freqframe']
398 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
399 'GEO','GALACTO','LGROUP','CMB']
400 if frame in valid:
401 inf = list(self._getcoordinfo())
402 inf[1] = frame
403 self._setcoordinfo(inf)
404 else:
405 print "Please specify a valid freq type. Valid types are:\n",valid
406
407 def get_unit(self):
408 """
409 Get the default unit set in this scantable
410 Parameters:
411 Returns:
412 A unit string
413 """
414 inf = self._getcoordinfo()
415 unit = inf[0]
416 if unit == '': unit = 'channel'
417 return unit
418
419 def get_abcissa(self, rowno=0):
420 """
421 Get the abcissa in the current coordinate setup for the currently
422 selected Beam/IF/Pol
423 Parameters:
424 rowno: an optional row number in the scantable. Default is the
425 first row, i.e. rowno=0
426 Returns:
427 The abcissa values and it's format string.
428 """
429 abc = self._getabcissa(rowno)
430 lbl = self._getabcissalabel(rowno)
431 return abc, lbl
432
433 def create_mask(self, *args, **kwargs):
434 """
435 Compute and return a mask based on [min,max] windows.
436 The specified windows are to be INCLUDED, when the mask is
437 applied.
438 Parameters:
439 [min,max],[min2,max2],...
440 Pairs of start/end points specifying the regions
441 to be masked
442 invert: optional argument. If specified as True,
443 return an inverted mask, i.e. the regions
444 specified are EXCLUDED
445 Example:
446 scan.set_unit('channel')
447
448 a)
449 msk = scan.set_mask([400,500],[800,900])
450 # masks everything outside 400 and 500
451 # and 800 and 900 in the unit 'channel'
452
453 b)
454 msk = scan.set_mask([400,500],[800,900], invert=True)
455 # masks the regions between 400 and 500
456 # and 800 and 900 in the unit 'channel'
457
458 """
459 u = self._getcoordinfo()[0]
460 if self._vb:
461 if u == "": u = "channel"
462 print "The current mask window unit is", u
463 n = self.nchan()
464 data = self._getabcissa()
465 msk = zeros(n)
466 for window in args:
467 if (len(window) != 2 or window[0] > window[1] ):
468 print "A window needs to be defined as [min,max]"
469 return
470 for i in range(n):
471 if data[i] >= window[0] and data[i] < window[1]:
472 msk[i] = 1
473 if kwargs.has_key('invert'):
474 if kwargs.get('invert'):
475 from numarray import logical_not
476 msk = logical_not(msk)
477 return msk
478
479 def set_restfreqs(self, freqs, unit='Hz'):
480 """
481 Set the restfrequency(s) for this scantable.
482 Parameters:
483 freqs: one or more frequencies
484 unit: optional 'unit', default 'Hz'
485 Example:
486 scan.set_restfreqs([1000000000.0])
487 """
488 if type(freqs) is float or int:
489 freqs = (freqs)
490 sdtable._setrestfreqs(self,freqs, unit)
491 return
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 flag_spectrum(self, thebeam, theif, thepol):
504 """
505 This flags a selected spectrum in the scan 'for good'.
506 USE WITH CARE - not reversible.
507 Use masks for non-permanent exclusion of channels.
508 Parameters:
509 thebeam,theif,thepol: all have to be explicitly
510 specified
511 Example:
512 scan.flag_spectrum(0,0,1)
513 flags the spectrum for Beam=0, IF=0, Pol=1
514 """
515 if (thebeam < self.nbeam() and theif < self.nif() and thepol < self.npol()):
516 stable.setbeam(thebeam)
517 stable.setif(theif)
518 stable.setpol(thepol)
519 stable._flag(self)
520 else:
521 print "Please specify a valid (Beam/IF/Pol)"
522 return
523
524 def plot(self, what='spectrum',col='Pol', panel=None):
525 """
526 Plot the spectra contained in the scan. Alternatively you can also
527 Plot Tsys vs Time
528 Parameters:
529 what: a choice of 'spectrum' (default) or 'tsys'
530 col: which out of Beams/IFs/Pols should be colour stacked
531 panel: set up multiple panels, currently not working.
532 """
533 print "Warning! Not fully functional. Use plotter.plot() instead"
534
535 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
536
537 validyax = ['spectrum','tsys']
538 from asap.asaplot import ASAPlot
539 if not self._p:
540 self._p = ASAPlot()
541 #print "Plotting not enabled"
542 #return
543 if self._p.is_dead:
544 del self._p
545 self._p = ASAPlot()
546 npan = 1
547 x = None
548 if what == 'tsys':
549 n = self.nrow()
550 if n < 2:
551 print "Only one integration. Can't plot."
552 return
553 self._p.hold()
554 self._p.clear()
555 if panel == 'Time':
556 npan = self.nrow()
557 self._p.set_panels(rows=npan)
558 xlab,ylab,tlab = None,None,None
559 self._vb = False
560 sel = self.get_cursor()
561 for i in range(npan):
562 if npan > 1:
563 self._p.subplot(i)
564 for j in range(validcol[col]):
565 x = None
566 y = None
567 m = None
568 tlab = self._getsourcename(i)
569 import re
570 tlab = re.sub('_S','',tlab)
571 if col == 'Beam':
572 self.setbeam(j)
573 elif col == 'IF':
574 self.setif(j)
575 elif col == 'Pol':
576 self.setpol(j)
577 if what == 'tsys':
578 x = range(self.nrow())
579 xlab = 'Time [pixel]'
580 m = list(ones(len(x)))
581 y = []
582 ylab = r'$T_{sys}$'
583 for k in range(len(x)):
584 y.append(self._gettsys(k))
585 else:
586 x,xlab = self.get_abcissa(i)
587 y = self._getspectrum(i)
588 ylab = r'Flux'
589 m = self._getmask(i)
590 llab = col+' '+str(j)
591 self._p.set_line(label=llab)
592 self._p.plot(x,y,m)
593 self._p.set_axes('xlabel',xlab)
594 self._p.set_axes('ylabel',ylab)
595 self._p.set_axes('title',tlab)
596 self._p.release()
597 self.set_cursor(sel[0],sel[1],sel[2])
598 self._vb = rcParams['verbose']
599 return
600
601 print out
602
603 def _print_values(self, dat, label='', timestamps=[]):
604 d = dat['data']
605 a = dat['axes']
606 shp = d.getshape()
607 out = ''
608 for i in range(shp[3]):
609 out += '%s [%s]:\n' % (a[3],timestamps[i])
610 t = d[:,:,:,i]
611 for j in range(shp[0]):
612 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
613 for k in range(shp[1]):
614 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
615 for l in range(shp[2]):
616 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
617 out += '= %3.3f\n' % (t[j,k,l])
618 out += "--------------------------------------------------\n"
619 print "--------------------------------------------------"
620 print " ", label
621 print "--------------------------------------------------"
622 print out
Note: See TracBrowser for help on using the repository browser.