source: trunk/python/scantable.py@ 351

Last change on this file since 351 was 346, checked in by kil064, 20 years ago

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