source: trunk/python/scantable.py@ 587

Last change on this file since 587 was 582, checked in by mar637, 20 years ago

fix for asap0003

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.6 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, average=None, 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 average: average all integrations withinb a scan on read.
24 The default (True) is taken from .asaprc.
25 unit: brightness unit; must be consistent with K or Jy.
26 Over-rides the default selected by the reader
27 (input rpfits/sdfits/ms) or replaces the value
28 in existing scantables
29 """
30 if average is None or type(average) is not bool:
31 average = rcParams['scantable.autoaverage']
32
33 varlist = vars()
34 self._vb = rcParams['verbose']
35 self._p = None
36
37 if isinstance(filename,sdtable):
38 sdtable.__init__(self, filename)
39 if unit is not None:
40 self.set_fluxunit(unit)
41 else:
42 import os.path
43 if not os.path.exists(filename):
44 print "File '%s' not found." % (filename)
45 return
46 filename = os.path.expandvars(filename)
47 if os.path.isdir(filename):
48 # crude check if asap table
49 if os.path.exists(filename+'/table.info'):
50 sdtable.__init__(self, filename)
51 if unit is not None:
52 self.set_fluxunit(unit)
53 else:
54 print "The given file '%s'is not a valid asap table." % (filename)
55 return
56 else:
57 from asap._asap import sdreader
58 ifSel = -1
59 beamSel = -1
60 r = sdreader(filename,ifSel,beamSel)
61 print 'Importing data...'
62 r._read([-1])
63 tbl = r._getdata()
64 if unit is not None:
65 tbl.set_fluxunit(unit)
66 if average:
67 from asap._asap import average as _av
68 print 'Auto averaging integrations...'
69 tbl2 = _av((tbl,),(),True,'none')
70 sdtable.__init__(self,tbl2)
71 del tbl2
72 else:
73 sdtable.__init__(self,tbl)
74 del r,tbl
75 self._add_history("scantable", varlist)
76
77 def save(self, name=None, format=None, stokes=False, overwrite=False):
78 """
79 Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
80 Image FITS or MS2 format.
81 Parameters:
82 name: the name of the outputfile. For format="FITS" this
83 is the directory file name into which all the files
84 will be written (default is 'asap_FITS'). For format
85 "ASCII" this is the root file name (data in 'name'.txt
86 and header in 'name'_header.txt)
87 format: an optional file format. Default is ASAP.
88 Allowed are - 'ASAP' (save as ASAP [aips++] Table),
89 'SDFITS' (save as SDFITS file)
90 'FITS' (saves each row as a FITS Image)
91 'ASCII' (saves as ascii text file)
92 'MS2' (saves as an aips++
93 MeasurementSet V2)
94 stokes: Convert to Stokes parameters (only available
95 currently with FITS and ASCII formats.
96 Default is False.
97 overwrite: If the file should be overwritten if it exists.
98 The default False is to return with warning
99 without writing the output. USE WITH CARE.
100 Example:
101 scan.save('myscan.asap')
102 scan.save('myscan.sdfits','SDFITS')
103 """
104 from os import path
105 if format is None: format = rcParams['scantable.save']
106 suffix = '.'+format.lower()
107 if name is None or name =="":
108 name = 'scantable'+suffix
109 print "No filename given. Using default name %s..." % name
110 name = path.expandvars(name)
111 if path.isfile(name) or path.isdir(name):
112 if not overwrite:
113 print "File %s already exists." % name
114 return
115 format2 = format.upper()
116 if format2 == 'ASAP':
117 self._save(name)
118 else:
119 from asap._asap import sdwriter as _sw
120 w = _sw(format2)
121 w.write(self, name, stokes)
122 return
123
124 def copy(self):
125 """
126 Return a copy of this scantable.
127 Parameters:
128 none
129 Example:
130 copiedscan = scan.copy()
131 """
132 sd = scantable(sdtable._copy(self))
133 return sd
134
135 def get_scan(self, scanid=None):
136 """
137 Return a specific scan (by scanno) or collection of scans (by
138 source name) in a new scantable.
139 Parameters:
140 scanid: a (list of) scanno or a source name, unix-style
141 patterns are accepted for source name matching, e.g.
142 '*_R' gets all 'ref scans
143 Example:
144 # get all scans containing the source '323p459'
145 newscan = scan.get_scan('323p459')
146 # get all 'off' scans
147 refscans = scan.get_scan('*_R')
148 # get a susbset of scans by scanno (as listed in scan.summary())
149 newscan = scan.get_scan([0,2,7,10])
150 """
151 if scanid is None:
152 print "Please specify a scan no or name to retrieve from the scantable"
153 try:
154 if type(scanid) is str:
155 s = sdtable._getsource(self,scanid)
156 return scantable(s)
157 elif type(scanid) is int:
158 s = sdtable._getscan(self,[scanid])
159 return scantable(s)
160 elif type(scanid) is list:
161 s = sdtable._getscan(self,scanid)
162 return scantable(s)
163 else:
164 print "Illegal scanid type, use 'int' or 'list' if ints."
165 except RuntimeError:
166 print "Couldn't find any match."
167
168 def __str__(self):
169 return sdtable._summary(self,True)
170
171 def summary(self,filename=None, verbose=None):
172 """
173 Print a summary of the contents of this scantable.
174 Parameters:
175 filename: the name of a file to write the putput to
176 Default - no file output
177 verbose: print extra info such as the frequency table
178 The default (False) is taken from .asaprc
179 """
180 info = sdtable._summary(self, verbose)
181 if verbose is None: verbose = rcParams['scantable.verbosesummary']
182 if filename is not None:
183 if filename is "":
184 filename = 'scantable_summary.txt'
185 from os.path import expandvars, isdir
186 filename = expandvars(filename)
187 if not isdir(filename):
188 data = open(filename, 'w')
189 data.write(info)
190 data.close()
191 else:
192 print "Illegal file name '%s'." % (filename)
193 print info
194
195 def set_cursor(self, beam=0, IF=0, pol=0):
196 """
197 Set the spectrum for individual operations.
198 Parameters:
199 beam, IF, pol: a number
200 Example:
201 scan.set_cursor(0,0,1)
202 pol1sig = scan.stats(all=False) # returns std dev for beam=0
203 # if=0, pol=1
204 """
205 varlist = vars()
206 self.setbeam(beam)
207 self.setpol(pol)
208 self.setif(IF)
209 self._add_history("set_cursor",varlist)
210 return
211
212 def get_cursor(self):
213 """
214 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
215 Parameters:
216 none
217 Returns:
218 a list of values (currentBeam,currentIF,currentPol)
219 Example:
220 none
221 """
222 i = self.getbeam()
223 j = self.getif()
224 k = self.getpol()
225 if self._vb:
226 print "--------------------------------------------------"
227 print " Cursor position"
228 print "--------------------------------------------------"
229 out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
230 print out
231 return i,j,k
232
233 def stats(self, stat='stddev', mask=None, allaxes=None):
234 """
235 Determine the specified statistic of the current beam/if/pol
236 Takes a 'mask' as an optional parameter to specify which
237 channels should be excluded.
238 Parameters:
239 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
240 'var', 'stddev', 'avdev', 'rms', 'median'
241 mask: an optional mask specifying where the statistic
242 should be determined.
243 allaxes: if True apply to all spectra. Otherwise
244 apply only to the selected (beam/pol/if)spectra only.
245 The default is taken from .asaprc (True if none)
246 Example:
247 scan.set_unit('channel')
248 msk = scan.create_mask([100,200],[500,600])
249 scan.stats(stat='mean', mask=m)
250 """
251 if allaxes is None: allaxes = rcParams['scantable.allaxes']
252 from asap._asap import stats as _stats
253 from numarray import array,zeros,Float
254 if mask == None:
255 mask = ones(self.nchan())
256 axes = ['Beam','IF','Pol','Time']
257
258 beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
259 if allaxes:
260 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
261 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
262 arr = array(zeros(n),shape=shp,type=Float)
263
264 for i in range(self.nbeam()):
265 self.setbeam(i)
266 for j in range(self.nif()):
267 self.setif(j)
268 for k in range(self.npol()):
269 self.setpol(k)
270 arr[i,j,k,:] = _stats(self,mask,stat,-1)
271 retval = {'axes': axes, 'data': arr, 'cursor':None}
272 tm = [self._gettime(val) for val in range(self.nrow())]
273 if self._vb:
274 self._print_values(retval,stat,tm)
275 self.setbeam(beamSel)
276 self.setif(IFSel)
277 self.setpol(polSel)
278 return retval
279
280 else:
281 statval = _stats(self,mask,stat,-1)
282 out = ''
283 for l in range(self.nrow()):
284 tm = self._gettime(l)
285 out += 'Time[%s]:\n' % (tm)
286 if self.nbeam() > 1: out += ' Beam[%d] ' % (beamSel)
287 if self.nif() > 1: out += ' IF[%d] ' % (IFSel)
288 if self.npol() > 1: out += ' Pol[%d] ' % (polSel)
289 out += '= %3.3f\n' % (statval[l])
290 out += "--------------------------------------------------\n"
291
292 if self._vb:
293 print "--------------------------------------------------"
294 print " ",stat
295 print "--------------------------------------------------"
296 print out
297 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
298 return retval
299
300 def stddev(self,mask=None, allaxes=None):
301 """
302 Determine the standard deviation of the current beam/if/pol
303 Takes a 'mask' as an optional parameter to specify which
304 channels should be excluded.
305 Parameters:
306 mask: an optional mask specifying where the standard
307 deviation should be determined.
308 allaxes: optional flag to show all or a cursor selected
309 spectrum of Beam/IF/Pol. Default is all or taken
310 from .asaprc
311
312 Example:
313 scan.set_unit('channel')
314 msk = scan.create_mask([100,200],[500,600])
315 scan.stddev(mask=m)
316 """
317 if allaxes is None: allaxes = rcParams['scantable.allaxes']
318 return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
319
320 def get_tsys(self, allaxes=None):
321 """
322 Return the System temperatures.
323 Parameters:
324 allaxes: if True apply to all spectra. Otherwise
325 apply only to the selected (beam/pol/if)spectra only.
326 The default is taken from .asaprc (True if none)
327 Returns:
328 a list of Tsys values.
329 """
330 if allaxes is None: allaxes = rcParams['scantable.allaxes']
331 from numarray import array,zeros,Float
332 axes = ['Beam','IF','Pol','Time']
333
334 if allaxes:
335 n = self.nbeam()*self.nif()*self.npol()*self.nrow()
336 shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
337 arr = array(zeros(n),shape=shp,type=Float)
338
339 for i in range(self.nbeam()):
340 self.setbeam(i)
341 for j in range(self.nif()):
342 self.setif(j)
343 for k in range(self.npol()):
344 self.setpol(k)
345 arr[i,j,k,:] = self._gettsys()
346 retval = {'axes': axes, 'data': arr, 'cursor':None}
347 tm = [self._gettime(val) for val in range(self.nrow())]
348 if self._vb:
349 self._print_values(retval,'Tsys',tm)
350 return retval
351
352 else:
353 i,j,k = (self.getbeam(),self.getif(),self.getpol())
354 statval = self._gettsys()
355 out = ''
356 for l in range(self.nrow()):
357 tm = self._gettime(l)
358 out += 'Time[%s]:\n' % (tm)
359 if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
360 if self.nif() > 1: out += ' IF[%d] ' % (j)
361 if self.npol() > 1: out += ' Pol[%d] ' % (k)
362 out += '= %3.3f\n' % (statval[l])
363 out += "--------------------------------------------------\n"
364
365 if self._vb:
366 print "--------------------------------------------------"
367 print " TSys"
368 print "--------------------------------------------------"
369 print out
370 retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
371 return retval
372
373 def get_time(self, row=-1):
374 """
375 Get a list of time stamps for the observations.
376 Return a string for each integration in the scantable.
377 Parameters:
378 row: row no of integration. Default -1 return all rows
379 Example:
380 none
381 """
382 out = []
383 if row == -1:
384 for i in range(self.nrow()):
385 out.append(self._gettime(i))
386 return out
387 else:
388 if row < self.nrow():
389 return self._gettime(row)
390
391 def set_unit(self, unit='channel'):
392 """
393 Set the unit for all following operations on this scantable
394 Parameters:
395 unit: optional unit, default is 'channel'
396 one of '*Hz','km/s','channel', ''
397 """
398 varlist = vars()
399 if unit in ['','pixel', 'channel']:
400 unit = ''
401 inf = list(self._getcoordinfo())
402 inf[0] = unit
403 self._setcoordinfo(inf)
404 if self._p: self.plot()
405 self._add_history("set_unit",varlist)
406
407 def set_instrument(self, instr):
408 """
409 Set the instrument for subsequent processing
410 Parameters:
411 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
412 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
413 """
414 self._setInstrument(instr)
415 self._add_history("set_instument",vars())
416
417 def set_doppler(self, doppler='RADIO'):
418 """
419 Set the doppler for all following operations on this scantable.
420 Parameters:
421 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
422 """
423 varlist = vars()
424 inf = list(self._getcoordinfo())
425 inf[2] = doppler
426 self._setcoordinfo(inf)
427 if self._p: self.plot()
428 self._add_history("set_doppler",vars())
429
430 def set_freqframe(self, frame=None):
431 """
432 Set the frame type of the Spectral Axis.
433 Parameters:
434 frame: an optional frame type, default 'LSRK'.
435 Examples:
436 scan.set_freqframe('BARY')
437 """
438 if frame is None: frame = rcParams['scantable.freqframe']
439 varlist = vars()
440 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
441 'GEO','GALACTO','LGROUP','CMB']
442 if frame in valid:
443 inf = list(self._getcoordinfo())
444 inf[1] = frame
445 self._setcoordinfo(inf)
446 self._add_history("set_freqframe",varlist)
447 else:
448 print "Please specify a valid freq type. Valid types are:\n",valid
449
450 def get_unit(self):
451 """
452 Get the default unit set in this scantable
453 Parameters:
454 Returns:
455 A unit string
456 """
457 inf = self._getcoordinfo()
458 unit = inf[0]
459 if unit == '': unit = 'channel'
460 return unit
461
462 def get_abcissa(self, rowno=0):
463 """
464 Get the abcissa in the current coordinate setup for the currently
465 selected Beam/IF/Pol
466 Parameters:
467 rowno: an optional row number in the scantable. Default is the
468 first row, i.e. rowno=0
469 Returns:
470 The abcissa values and it's format string (as a dictionary)
471 """
472 abc = self._getabcissa(rowno)
473 lbl = self._getabcissalabel(rowno)
474 return abc, lbl
475 #return {'abcissa':abc,'label':lbl}
476
477 def create_mask(self, *args, **kwargs):
478 """
479 Compute and return a mask based on [min,max] windows.
480 The specified windows are to be INCLUDED, when the mask is
481 applied.
482 Parameters:
483 [min,max],[min2,max2],...
484 Pairs of start/end points specifying the regions
485 to be masked
486 invert: optional argument. If specified as True,
487 return an inverted mask, i.e. the regions
488 specified are EXCLUDED
489 row: create the mask using the specified row for
490 unit conversions, default is row=0
491 only necessary if frequency varies over rows.
492 Example:
493 scan.set_unit('channel')
494
495 a)
496 msk = scan.set_mask([400,500],[800,900])
497 # masks everything outside 400 and 500
498 # and 800 and 900 in the unit 'channel'
499
500 b)
501 msk = scan.set_mask([400,500],[800,900], invert=True)
502 # masks the regions between 400 and 500
503 # and 800 and 900 in the unit 'channel'
504
505 """
506 row = 0
507 if kwargs.has_key("row"):
508 row = kwargs.get("row")
509 data = self._getabcissa(row)
510 u = self._getcoordinfo()[0]
511 if self._vb:
512 if u == "": u = "channel"
513 print "The current mask window unit is", u
514 n = self.nchan()
515 msk = zeros(n)
516 for window in args:
517 if (len(window) != 2 or window[0] > window[1] ):
518 print "A window needs to be defined as [min,max]"
519 return
520 for i in range(n):
521 if data[i] >= window[0] and data[i] < window[1]:
522 msk[i] = 1
523 if kwargs.has_key('invert'):
524 if kwargs.get('invert'):
525 from numarray import logical_not
526 msk = logical_not(msk)
527 return msk
528
529 def get_restfreqs(self):
530 """
531 Get the restfrequency(s) stored in this scantable.
532 The return value(s) are always of unit 'Hz'
533 Parameters:
534 none
535 Returns:
536 a list of doubles
537 """
538 return list(self._getrestfreqs())
539
540 def lines(self):
541 """
542 Print the list of known spectral lines
543 """
544 sdtable._lines(self)
545
546 def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
547 theif=None):
548 """
549 Select the restfrequency for the specified source and IF OR
550 replace for all IFs. If the 'freqs' argument holds a scalar,
551 then that rest frequency will be applied to the selected
552 data (and added to the list of available rest frequencies).
553 In this way, you can set a rest frequency for each
554 source and IF combination. If the 'freqs' argument holds
555 a vector, then it MUST be of length the number of IFs
556 (and the available restfrequencies will be replaced by
557 this vector). In this case, *all* data ('source' and
558 'theif' are ignored) have the restfrequency set per IF according
559 to the corresponding value you give in the 'freqs' vector.
560 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
561 IF 1 gets restfreq 2e9.
562
563 You can also specify the frequencies via known line names
564 in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
565 takes precedence. See the list of known names via function
566 scantable.lines()
567 Parameters:
568 freqs: list of rest frequencies
569 unit: unit for rest frequency (default 'Hz')
570 lines: list of known spectral lines names (alternative to freqs).
571 See possible list via scantable.lines()
572 source: Source name (blank means all)
573 theif: IF (-1 means all)
574 Example:
575 scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
576 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
577 """
578 varlist = vars()
579 if source is None:
580 source = ""
581 if theif is None:
582 theif = -1
583 t = type(freqs)
584 if t is int or t is float:
585 freqs = [freqs]
586 if freqs is None:
587 freqs = []
588 t = type(lines)
589 if t is str:
590 lines = [lines]
591 if lines is None:
592 lines = []
593 sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
594 self._add_history("set_restfreqs", varlist)
595
596
597
598 def flag_spectrum(self, thebeam, theif, thepol):
599 """
600 This flags a selected spectrum in the scan 'for good'.
601 USE WITH CARE - not reversible.
602 Use masks for non-permanent exclusion of channels.
603 Parameters:
604 thebeam,theif,thepol: all have to be explicitly
605 specified
606 Example:
607 scan.flag_spectrum(0,0,1)
608 flags the spectrum for Beam=0, IF=0, Pol=1
609 """
610 if (thebeam < self.nbeam() and
611 theif < self.nif() and
612 thepol < self.npol()):
613 sdtable.setbeam(self, thebeam)
614 sdtable.setif(self, theif)
615 sdtable.setpol(self, thepol)
616 sdtable._flag(self)
617 self._add_history("flag_spectrum", vars())
618 else:
619 print "Please specify a valid (Beam/IF/Pol)"
620 return
621
622 def plot(self, what='spectrum',col='Pol', panel=None):
623 """
624 Plot the spectra contained in the scan. Alternatively you can also
625 Plot Tsys vs Time
626 Parameters:
627 what: a choice of 'spectrum' (default) or 'tsys'
628 col: which out of Beams/IFs/Pols should be colour stacked
629 panel: set up multiple panels, currently not working.
630 """
631 print "Warning! Not fully functional. Use plotter.plot() instead"
632
633 validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
634
635 validyax = ['spectrum','tsys']
636 from asap.asaplot import ASAPlot
637 if not self._p:
638 self._p = ASAPlot()
639 #print "Plotting not enabled"
640 #return
641 if self._p.is_dead:
642 del self._p
643 self._p = ASAPlot()
644 npan = 1
645 x = None
646 if what == 'tsys':
647 n = self.nrow()
648 if n < 2:
649 print "Only one integration. Can't plot."
650 return
651 self._p.hold()
652 self._p.clear()
653 if panel == 'Time':
654 npan = self.nrow()
655 self._p.set_panels(rows=npan)
656 xlab,ylab,tlab = None,None,None
657 self._vb = False
658 sel = self.get_cursor()
659 for i in range(npan):
660 if npan > 1:
661 self._p.subplot(i)
662 for j in range(validcol[col]):
663 x = None
664 y = None
665 m = None
666 tlab = self._getsourcename(i)
667 import re
668 tlab = re.sub('_S','',tlab)
669 if col == 'Beam':
670 self.setbeam(j)
671 elif col == 'IF':
672 self.setif(j)
673 elif col == 'Pol':
674 self.setpol(j)
675 if what == 'tsys':
676 x = range(self.nrow())
677 xlab = 'Time [pixel]'
678 m = list(ones(len(x)))
679 y = []
680 ylab = r'$T_{sys}$'
681 for k in range(len(x)):
682 y.append(self._gettsys(k))
683 else:
684 x,xlab = self.get_abcissa(i)
685 y = self._getspectrum(i)
686 ylab = r'Flux'
687 m = self._getmask(i)
688 llab = col+' '+str(j)
689 self._p.set_line(label=llab)
690 self._p.plot(x,y,m)
691 self._p.set_axes('xlabel',xlab)
692 self._p.set_axes('ylabel',ylab)
693 self._p.set_axes('title',tlab)
694 self._p.release()
695 self.set_cursor(sel[0],sel[1],sel[2])
696 self._vb = rcParams['verbose']
697 return
698
699 print out
700
701 def _print_values(self, dat, label='', timestamps=[]):
702 d = dat['data']
703 a = dat['axes']
704 shp = d.getshape()
705 out = ''
706 for i in range(shp[3]):
707 out += '%s [%s]:\n' % (a[3],timestamps[i])
708 t = d[:,:,:,i]
709 for j in range(shp[0]):
710 if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
711 for k in range(shp[1]):
712 if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
713 for l in range(shp[2]):
714 if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
715 out += '= %3.3f\n' % (t[j,k,l])
716 out += "-"*80
717 out += "\n"
718 print "-"*80
719 print " ", label
720 print "-"*80
721 print out
722
723 def history(self):
724 hist = list(self._gethistory())
725 print "-"*80
726 for h in hist:
727 if h.startswith("---"):
728 print h
729 else:
730 items = h.split("##")
731 date = items[0]
732 func = items[1]
733 items = items[2:]
734 print date
735 print "Function: %s\n Parameters:" % (func)
736 for i in items:
737 s = i.split("=")
738 print " %s = %s" % (s[0],s[1])
739 print "-"*80
740 return
741
742 #
743 # Maths business
744 #
745
746 def average_time(self, mask=None, scanav=False, weight='tint'):
747 """
748 Return the (time) average of a scan, or apply it 'insitu'.
749 Note:
750 in channels only
751 The cursor of the output scan is set to 0.
752 Parameters:
753 one scan or comma separated scans
754 mask: an optional mask (only used for 'var' and 'tsys'
755 weighting)
756 scanav: True averages each scan separately
757 False (default) averages all scans together,
758 weight: Weighting scheme. 'none', 'var' (1/var(spec)
759 weighted), 'tsys' (1/Tsys**2 weighted), 'tint'
760 (integration time weighted) or 'tintsys' (Tint/Tsys**2).
761 The default is 'tint'
762 Example:
763 # time average the scantable without using a mask
764 newscan = scan.average_time()
765 """
766 varlist = vars()
767 if weight is None: weight = 'tint'
768 if mask is None: mask = ()
769 from asap._asap import average as _av
770 s = scantable(_av((self,), mask, scanav, weight))
771 s._add_history("average_time",varlist)
772 return s
773
774 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
775 allaxes=None):
776 """
777 Return a scan where all spectra are converted to either
778 Jansky or Kelvin depending upon the flux units of the scan table.
779 By default the function tries to look the values up internally.
780 If it can't find them (or if you want to over-ride), you must
781 specify EITHER jyperk OR eta (and D which it will try to look up
782 also if you don't set it). jyperk takes precedence if you set both.
783 Parameters:
784 jyperk: the Jy / K conversion factor
785 eta: the aperture efficiency
786 d: the geomtric diameter (metres)
787 insitu: if False a new scantable is returned.
788 Otherwise, the scaling is done in-situ
789 The default is taken from .asaprc (False)
790 allaxes: if True apply to all spectra. Otherwise
791 apply only to the selected (beam/pol/if)spectra only
792 The default is taken from .asaprc (True if none)
793 """
794 if allaxes is None: allaxes = rcParams['scantable.allaxes']
795 if insitu is None: insitu = rcParams['insitu']
796 varlist = vars()
797 if jyperk is None: jyperk = -1.0
798 if d is None: d = -1.0
799 if eta is None: eta = -1.0
800 if not insitu:
801 from asap._asap import convertflux as _convert
802 s = scantable(_convert(self, d, eta, jyperk, allaxes))
803 s._add_history("convert_flux", varlist)
804 return s
805 else:
806 from asap._asap import convertflux_insitu as _convert
807 _convert(self, d, eta, jyperk, allaxes)
808 self._add_history("convert_flux", varlist)
809 return
810
811 def gain_el(self, poly=None, filename="", method="linear",
812 insitu=None, allaxes=None):
813 """
814 Return a scan after applying a gain-elevation correction.
815 The correction can be made via either a polynomial or a
816 table-based interpolation (and extrapolation if necessary).
817 You specify polynomial coefficients, an ascii table or neither.
818 If you specify neither, then a polynomial correction will be made
819 with built in coefficients known for certain telescopes (an error
820 will occur if the instrument is not known).
821 The data and Tsys are *divided* by the scaling factors.
822 Parameters:
823 poly: Polynomial coefficients (default None) to compute a
824 gain-elevation correction as a function of
825 elevation (in degrees).
826 filename: The name of an ascii file holding correction factors.
827 The first row of the ascii file must give the column
828 names and these MUST include columns
829 "ELEVATION" (degrees) and "FACTOR" (multiply data
830 by this) somewhere.
831 The second row must give the data type of the
832 column. Use 'R' for Real and 'I' for Integer.
833 An example file would be
834 (actual factors are arbitrary) :
835
836 TIME ELEVATION FACTOR
837 R R R
838 0.1 0 0.8
839 0.2 20 0.85
840 0.3 40 0.9
841 0.4 60 0.85
842 0.5 80 0.8
843 0.6 90 0.75
844 method: Interpolation method when correcting from a table.
845 Values are "nearest", "linear" (default), "cubic"
846 and "spline"
847 insitu: if False a new scantable is returned.
848 Otherwise, the scaling is done in-situ
849 The default is taken from .asaprc (False)
850 allaxes: If True apply to all spectra. Otherwise
851 apply only to the selected (beam/pol/if) spectra only
852 The default is taken from .asaprc (True if none)
853 """
854
855 if allaxes is None: allaxes = rcParams['scantable.allaxes']
856 if insitu is None: insitu = rcParams['insitu']
857 varlist = vars()
858 if poly is None:
859 poly = ()
860 from os.path import expandvars
861 filename = expandvars(filename)
862 if not insitu:
863 from asap._asap import gainel as _gainEl
864 s = scantable(_gainEl(self, poly, filename, method, allaxes))
865 s._add_history("gain_el", varlist)
866 return s
867 else:
868 from asap._asap import gainel_insitu as _gainEl
869 _gainEl(self, poly, filename, method, allaxes)
870 self._add_history("gain_el", varlist)
871 return
872
873 def freq_align(self, reftime=None, method='cubic', perif=False,
874 insitu=None):
875 """
876 Return a scan where all rows have been aligned in frequency/velocity.
877 The alignment frequency frame (e.g. LSRK) is that set by function
878 set_freqframe.
879 Parameters:
880 reftime: reference time to align at. By default, the time of
881 the first row of data is used.
882 method: Interpolation method for regridding the spectra.
883 Choose from "nearest", "linear", "cubic" (default)
884 and "spline"
885 perif: Generate aligners per freqID (no doppler tracking) or
886 per IF (scan-based doppler tracking)
887 insitu: if False a new scantable is returned.
888 Otherwise, the scaling is done in-situ
889 The default is taken from .asaprc (False)
890 """
891 if insitu is None: insitu = rcParams['insitu']
892 varlist = vars()
893 if reftime is None: reftime = ''
894 perfreqid = not perif
895 if not insitu:
896 from asap._asap import freq_align as _align
897 s = scantable(_align(self, reftime, method, perfreqid))
898 s._add_history("freq_align", varlist)
899 return s
900 else:
901 from asap._asap import freq_align_insitu as _align
902 _align(self, reftime, method, perfreqid)
903 self._add_history("freq_align", varlist)
904 return
905
906 def opacity(self, tau, insitu=None, allaxes=None):
907 """
908 Apply an opacity correction. The data
909 and Tsys are multiplied by the correction factor.
910 Parameters:
911 tau: Opacity from which the correction factor is
912 exp(tau*ZD)
913 where ZD is the zenith-distance
914 insitu: if False a new scantable is returned.
915 Otherwise, the scaling is done in-situ
916 The default is taken from .asaprc (False)
917 allaxes: if True apply to all spectra. Otherwise
918 apply only to the selected (beam/pol/if)spectra only
919 The default is taken from .asaprc (True if none)
920 """
921 if allaxes is None: allaxes = rcParams['scantable.allaxes']
922 if insitu is None: insitu = rcParams['insitu']
923 varlist = vars()
924 if not insitu:
925 from asap._asap import opacity as _opacity
926 s = scantable(_opacity(self, tau, allaxes))
927 s._add_history("opacity", varlist)
928 return s
929 else:
930 from asap._asap import opacity_insitu as _opacity
931 _opacity(self, tau, allaxes)
932 self._add_history("opacity", varlist)
933 return
934
935 def bin(self, width=5, insitu=None):
936 """
937 Return a scan where all spectra have been binned up.
938 width: The bin width (default=5) in pixels
939 insitu: if False a new scantable is returned.
940 Otherwise, the scaling is done in-situ
941 The default is taken from .asaprc (False)
942 """
943 if insitu is None: insitu = rcParams['insitu']
944 varlist = vars()
945 if not insitu:
946 from asap._asap import bin as _bin
947 s = scantable(_bin(self, width))
948 s._add_history("bin",varlist)
949 return s
950 else:
951 from asap._asap import bin_insitu as _bin
952 _bin(self, width)
953 self._add_history("bin",varlist)
954 return
955
956
957 def resample(self, width=5, method='cubic', insitu=None):
958 """
959 Return a scan where all spectra have been binned up
960 width: The bin width (default=5) in pixels
961 method: Interpolation method when correcting from a table.
962 Values are "nearest", "linear", "cubic" (default)
963 and "spline"
964 insitu: if False a new scantable is returned.
965 Otherwise, the scaling is done in-situ
966 The default is taken from .asaprc (False)
967 """
968 if insitu is None: insitu = rcParams['insitu']
969 varlist = vars()
970 if not insitu:
971 from asap._asap import resample as _resample
972 s = scantable(_resample(self, method, width))
973 s._add_history("resample",varlist)
974 return s
975 else:
976 from asap._asap import resample_insitu as _resample
977 _resample(self, method, width)
978 self._add_history("resample",varlist)
979 return
980
981 def average_pol(self, mask=None, weight='none', insitu=None):
982 """
983 Average the Polarisations together.
984 The polarisation cursor of the output scan is set to 0
985 Parameters:
986 scan: The scantable
987 mask: An optional mask defining the region, where the
988 averaging will be applied. The output will have all
989 specified points masked.
990 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
991 weighted), or 'tsys' (1/Tsys**2 weighted)
992 insitu: if False a new scantable is returned.
993 Otherwise, the scaling is done in-situ
994 The default is taken from .asaprc (False)
995 Example:
996 polav = average_pols(myscan)
997 """
998 if insitu is None: insitu = rcParams['insitu']
999 varlist = vars()
1000
1001 if mask is None:
1002 mask = ()
1003 if not insitu:
1004 from asap._asap import averagepol as _avpol
1005 s = scantable(_avpol(self, mask, weight))
1006 s._add_history("average_pol",varlist)
1007 return s
1008 else:
1009 from asap._asap import averagepol_insitu as _avpol
1010 _avpol(self, mask, weight)
1011 self._add_history("average_pol",varlist)
1012 return
1013
1014 def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
1015 """
1016 Smooth the spectrum by the specified kernel (conserving flux).
1017 Parameters:
1018 scan: The input scan
1019 kernel: The type of smoothing kernel. Select from
1020 'hanning' (default), 'gaussian' and 'boxcar'.
1021 The first three characters are sufficient.
1022 width: The width of the kernel in pixels. For hanning this is
1023 ignored otherwise it defauls to 5 pixels.
1024 For 'gaussian' it is the Full Width Half
1025 Maximum. For 'boxcar' it is the full width.
1026 insitu: if False a new scantable is returned.
1027 Otherwise, the scaling is done in-situ
1028 The default is taken from .asaprc (False)
1029 allaxes: If True (default) apply to all spectra. Otherwise
1030 apply only to the selected (beam/pol/if)spectra only
1031 The default is taken from .asaprc (True if none)
1032 Example:
1033 none
1034 """
1035 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1036 if insitu is None: insitu = rcParams['insitu']
1037 varlist = vars()
1038 if not insitu:
1039 from asap._asap import smooth as _smooth
1040 s = scantable(_smooth(self,kernel,width,allaxes))
1041 s._add_history("smooth", varlist)
1042 return s
1043 else:
1044 from asap._asap import smooth_insitu as _smooth
1045 _smooth(self,kernel,width,allaxes)
1046 self._add_history("smooth", varlist)
1047 return
1048
1049 def poly_baseline(self, mask=None, order=0, insitu=None):
1050 """
1051 Return a scan which has been baselined (all rows) by a polynomial.
1052 Parameters:
1053 scan: a scantable
1054 mask: an optional mask
1055 order: the order of the polynomial (default is 0)
1056 insitu: if False a new scantable is returned.
1057 Otherwise, the scaling is done in-situ
1058 The default is taken from .asaprc (False)
1059 Example:
1060 # return a scan baselined by a third order polynomial,
1061 # not using a mask
1062 bscan = scan.poly_baseline(order=3)
1063 """
1064 if insitu is None: insitu = rcParams['insitu']
1065 varlist = vars()
1066 if mask is None:
1067 from numarray import ones
1068 mask = list(ones(scan.nchan()))
1069 from asap.asapfitter import fitter
1070 f = fitter()
1071 f._verbose(True)
1072 f.set_scan(self, mask)
1073 f.set_function(poly=order)
1074 sf = f.auto_fit(insitu)
1075 if insitu:
1076 self._add_history("poly_baseline", varlist)
1077 return
1078 else:
1079 sf._add_history("poly_baseline", varlist)
1080 return sf
1081
1082 def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
1083 threshold=3,insitu=None):
1084 """
1085 Return a scan which has been baselined (all rows) by a polynomial.
1086 Spectral lines are detected first using linefinder and masked out
1087 to avoid them affecting the baseline solution.
1088
1089 Parameters:
1090 scan: a scantable
1091 mask: an optional mask retreived from scantable
1092 edge: an optional number of channel to drop at
1093 the edge of spectrum. If only one value is
1094 specified, the same number will be dropped from
1095 both sides of the spectrum. Default is to keep
1096 all channels
1097 order: the order of the polynomial (default is 0)
1098 threshold: the threshold used by line finder. It is better to
1099 keep it large as only strong lines affect the
1100 baseline solution.
1101 insitu: if False a new scantable is returned.
1102 Otherwise, the scaling is done in-situ
1103 The default is taken from .asaprc (False)
1104
1105 Example:
1106 scan2=scan.auto_poly_baseline(order=7)
1107 """
1108 if insitu is None: insitu = rcParams['insitu']
1109 varlist = vars()
1110 from asap.asapfitter import fitter
1111 from asap.asaplinefind import linefinder
1112 from asap import _is_sequence_or_number as _is_valid
1113
1114 if not _is_valid(edge, int):
1115 raise RuntimeError, "Parameter 'edge' has to be an integer or a \
1116 pair of integers specified as a tuple"
1117
1118 # setup fitter
1119 f = fitter()
1120 f._verbose(True)
1121 f.set_function(poly=order)
1122
1123 # setup line finder
1124 fl=linefinder()
1125 fl.set_options(threshold=threshold)
1126
1127 if not insitu:
1128 workscan=self.copy()
1129 else:
1130 workscan=self
1131
1132 vb=workscan._vb
1133 # remember the verbose parameter and selection
1134 workscan._vb=False
1135 sel=workscan.get_cursor()
1136 rows=range(workscan.nrow())
1137 for i in range(workscan.nbeam()):
1138 workscan.setbeam(i)
1139 for j in range(workscan.nif()):
1140 workscan.setif(j)
1141 for k in range(workscan.npol()):
1142 workscan.setpol(k)
1143 if f._vb:
1144 print "Processing:"
1145 print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
1146 for iRow in rows:
1147 fl.set_scan(workscan,mask,edge)
1148 fl.find_lines(iRow)
1149 f.set_scan(workscan, fl.get_mask())
1150 f.x=workscan._getabcissa(iRow)
1151 f.y=workscan._getspectrum(iRow)
1152 f.data=None
1153 f.fit()
1154 x=f.get_parameters()
1155 workscan._setspectrum(f.fitter.getresidual(),iRow)
1156 workscan.set_cursor(sel[0],sel[1],sel[2])
1157 workscan._vb = vb
1158 if not insitu:
1159 return workscan
1160
1161 def rotate_linpolphase(self, angle, allaxes=None):
1162 """
1163 Rotate the phase of the complex polarization O=Q+iU correlation.
1164 This is always done in situ in the raw data. So if you call this
1165 function more than once then each call rotates the phase further.
1166 Parameters:
1167 angle: The angle (degrees) to rotate (add) by.
1168 allaxes: If True apply to all spectra. Otherwise
1169 apply only to the selected (beam/pol/if)spectra only.
1170 The default is taken from .asaprc (True if none)
1171 Examples:
1172 scan.rotate_linpolphase(2.3)
1173 """
1174 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1175 varlist = vars()
1176 from asap._asap import _rotate_linpolphase as _rotate
1177 _rotate(self, angle, allaxes)
1178 self._add_history("rotate_linpolphase", varlist)
1179 return
1180
1181
1182 def rotate_xyphase(self, angle, allaxes=None):
1183 """
1184 Rotate the phase of the XY correlation. This is always done in situ
1185 in the data. So if you call this function more than once
1186 then each call rotates the phase further.
1187 Parameters:
1188 angle: The angle (degrees) to rotate (add) by.
1189 allaxes: If True apply to all spectra. Otherwise
1190 apply only to the selected (beam/pol/if)spectra only.
1191 The default is taken from .asaprc (True if none)
1192 Examples:
1193 scan.rotate_xyphase(2.3)
1194 """
1195 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1196 varlist = vars()
1197 from asap._asap import rotate_xyphase as _rotate_xyphase
1198 _rotate_xyphase(self, angle, allaxes)
1199 self._add_history("rotate_xyphase", varlist)
1200 return
1201
1202
1203 def add(self, offset, insitu=None, allaxes=None):
1204 """
1205 Return a scan where all spectra have the offset added
1206 Parameters:
1207 offset: the offset
1208 insitu: if False a new scantable is returned.
1209 Otherwise, the scaling is done in-situ
1210 The default is taken from .asaprc (False)
1211 allaxes: if True apply to all spectra. Otherwise
1212 apply only to the selected (beam/pol/if)spectra only
1213 The default is taken from .asaprc (True if none)
1214 """
1215 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1216 if insitu is None: insitu = rcParams['insitu']
1217 varlist = vars()
1218 if not insitu:
1219 from asap._asap import add as _add
1220 s = scantable(_add(self, offset, allaxes))
1221 s._add_history("add",varlist)
1222 return s
1223 else:
1224 from asap._asap import add_insitu as _add
1225 _add(self, offset, allaxes)
1226 self._add_history("add",varlist)
1227 return
1228
1229 def scale(self, factor, insitu=None, allaxes=None, tsys=True):
1230 """
1231 Return a scan where all spectra are scaled by the give 'factor'
1232 Parameters:
1233 factor: the scaling factor
1234 insitu: if False a new scantable is returned.
1235 Otherwise, the scaling is done in-situ
1236 The default is taken from .asaprc (False)
1237 allaxes: if True apply to all spectra. Otherwise
1238 apply only to the selected (beam/pol/if)spectra only.
1239 The default is taken from .asaprc (True if none)
1240 tsys: if True (default) then apply the operation to Tsys
1241 as well as the data
1242 """
1243 if allaxes is None: allaxes = rcParams['scantable.allaxes']
1244 if insitu is None: insitu = rcParams['insitu']
1245 varlist = vars()
1246 if not insitu:
1247 from asap._asap import scale as _scale
1248 s = scantable(_scale(self, factor, allaxes, tsys))
1249 s._add_history("scale",varlist)
1250 return s
1251 else:
1252 from asap._asap import scale_insitu as _scale
1253 _scale(self, factor, allaxes, tsys)
1254 self._add_history("scale",varlist)
1255 return
1256
1257
1258 def quotient(self, other, isreference=True, preserve=True):
1259 """
1260 Return the quotient of a 'source' (on) scan and a 'reference' (off)
1261 scan.
1262 The reference can have just one row, even if the signal has many.
1263 Otherwise they must have the same number of rows.
1264 The cursor of the output scan is set to 0
1265 Parameters:
1266 other: the 'other' scan
1267 isreference: if the 'other' scan is the reference (default)
1268 or source
1269 preserve: you can preserve (default) the continuum or
1270 remove it. The equations used are
1271 preserve: Output = Toff * (on/off) - Toff
1272 remove: Output = Tref * (on/off) - Ton
1273 Example:
1274 # src is a scantable for an 'on' scan, ref for an 'off' scantable
1275 q1 = src.quotient(ref)
1276 q2 = ref.quotient(src, isreference=False)
1277 # gives the same result
1278 """
1279 from asap._asap import quotient as _quot
1280 if isreference:
1281 return scantable(_quot(self, other, preserve))
1282 else:
1283 return scantable(_quot(other, self, preserve))
1284
1285 def __add__(self, other):
1286 varlist = vars()
1287 s = None
1288 if isinstance(other, scantable):
1289 from asap._asap import b_operate as _bop
1290 s = scantable(_bop(self, other, 'add', True))
1291 elif isinstance(other, float):
1292 from asap._asap import add as _add
1293 s = scantable(_add(self, other, True))
1294 else:
1295 print "Other input is not a scantable or float value"
1296 return
1297 s._add_history("operator +", varlist)
1298 return s
1299
1300
1301 def __sub__(self, other):
1302 """
1303 implicit on all axes and on Tsys
1304 """
1305 varlist = vars()
1306 s = None
1307 if isinstance(other, scantable):
1308 from asap._asap import b_operate as _bop
1309 s = scantable(_bop(self, other, 'sub', True))
1310 elif isinstance(other, float):
1311 from asap._asap import add as _add
1312 s = scantable(_add(self, -other, True))
1313 else:
1314 print "Other input is not a scantable or float value"
1315 return
1316 s._add_history("operator -", varlist)
1317 return s
1318
1319 def __mul__(self, other):
1320 """
1321 implicit on all axes and on Tsys
1322 """
1323 varlist = vars()
1324 s = None
1325 if isinstance(other, scantable):
1326 from asap._asap import b_operate as _bop
1327 s = scantable(_bop(self, other, 'mul', True))
1328 elif isinstance(other, float):
1329 if other == 0.0:
1330 print "Multiplying by zero is not recommended."
1331 return
1332 from asap._asap import scale as _sca
1333 s = scantable(_sca(self, other, True))
1334 else:
1335 print "Other input is not a scantable or float value"
1336 return
1337 s._add_history("operator *", varlist)
1338 return s
1339
1340
1341 def __div__(self, other):
1342 """
1343 implicit on all axes and on Tsys
1344 """
1345 varlist = vars()
1346 s = None
1347 if isinstance(other, scantable):
1348 from asap._asap import b_operate as _bop
1349 s = scantable(_bop(self, other, 'div', True))
1350 elif isinstance(other, float):
1351 if other == 0.0:
1352 print "Dividing by zero is not recommended"
1353 return
1354 from asap._asap import scale as _sca
1355 s = scantable(_sca(self, 1.0/other, True))
1356 else:
1357 print "Other input is not a scantable or float value"
1358 return
1359 s._add_history("operator /", varlist)
1360 return s
1361
1362 def get_fit(self, row=0):
1363 """
1364 Print or return the stored fits for a row in the scantable
1365 Parameters:
1366 row: the row which the fit has been applied to.
1367 """
1368 if row > self.nrow():
1369 return
1370 from asap import asapfit
1371 fit = asapfit(self._getfit(row))
1372 if self._vb:
1373 print fit
1374 return
1375 else:
1376 return fit.as_dict()
1377
1378 def _add_history(self, funcname, parameters):
1379 # create date
1380 sep = "##"
1381 from datetime import datetime
1382 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1383 hist = dstr+sep
1384 hist += funcname+sep#cdate+sep
1385 if parameters.has_key('self'): del parameters['self']
1386 for k,v in parameters.iteritems():
1387 if type(v) is dict:
1388 for k2,v2 in v.iteritems():
1389 hist += k2
1390 hist += "="
1391 if isinstance(v2,scantable):
1392 hist += 'scantable'
1393 elif k2 == 'mask':
1394 if isinstance(v2,list) or isinstance(v2,tuple):
1395 hist += str(self._zip_mask(v2))
1396 else:
1397 hist += str(v2)
1398 else:
1399 hist += str(v2)
1400 else:
1401 hist += k
1402 hist += "="
1403 if isinstance(v,scantable):
1404 hist += 'scantable'
1405 elif k == 'mask':
1406 if isinstance(v,list) or isinstance(v,tuple):
1407 hist += str(self._zip_mask(v))
1408 else:
1409 hist += str(v)
1410 else:
1411 hist += str(v)
1412 hist += sep
1413 hist = hist[:-2] # remove trailing '##'
1414 self._addhistory(hist)
1415
1416
1417 def _zip_mask(self, mask):
1418 mask = list(mask)
1419 i = 0
1420 segments = []
1421 while mask[i:].count(1):
1422 i += mask[i:].index(1)
1423 if mask[i:].count(0):
1424 j = i + mask[i:].index(0)
1425 else:
1426 j = len(mask)
1427 segments.append([i,j])
1428 i = j
1429 return segments
Note: See TracBrowser for help on using the repository browser.