source: trunk/python/scantable.py@ 412

Last change on this file since 412 was 411, checked in by mar637, 20 years ago

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